1(in-package "FRICAS-LISP") 2 3(locally 4 (declare (optimize (speed 3) (safety 0))) 5 6#+(and :openmcl (or :X8632-TARGET :X8664-TARGET :ARM-TARGET)) 7(macrolet ( 8 (bignum_subtag() 9 #+:ARM-TARGET 'ARM::SUBTAG-BIGNUM 10 #+:PPC32-TARGET 'PPC32::SUBTAG-BIGNUM 11 #+:PPC64-TARGET 'PPC64::SUBTAG-BIGNUM 12 #+:X8632-TARGET 'X8632::SUBTAG-BIGNUM 13 #+:X8664-TARGET 'X8664::SUBTAG-BIGNUM 14 ) 15 (digits_to_words(dl) 16 #+:32-BIT-TARGET dl 17 #+:64-BIT-TARGET `(ceiling ,dl 2)) 18 (words_to_digits(wl) 19 #+:32-BIT-TARGET wl 20 #+:64-BIT-TARGET `(ash ,wl 1)) 21 (words_to_bytes(wl) 22 #+:32-BIT-TARGET `(ash ,wl 2) 23 #+:64-BIT-TARGET `(ash ,wl 3)) 24 (is_plus(n nl) 25 #-:X8664-TARGET `(eql (the fixnum (ccl::%bignum-sign ,n)) 0) 26 #+:X8664-TARGET `(ccl::%bignum-0-or-plusp ,n ,nl)) 27 (alloc_bignum (rl res) 28 #+(and :64-BIT-TARGET :CCL-1.12) 29 `(ccl::%maybe-allocate-bignum ,rl ,res) 30 #-(and :64-BIT-TARGET :CCL-1.12) 31 `(ccl::%alloc-misc ,rl (bignum_subtag))) 32 (maybe_call_with_res_2 (f x y) 33 #+(and :64-BIT-TARGET :CCL-1.12) 34 `(,f ,x ,y res) 35 #-(and :64-BIT-TARGET :CCL-1.12) 36 `(,f ,x ,y)) 37 ) 38 39#+:X8632-TARGET 40(progn 41(CCL::defx8632lapfunction gmp-bignum-copy-from-lisp 42 ((x 4) #|(ra 0)|# (y arg_y) (l arg_z)) 43 (mark-as-imm temp1) 44 (movl (@ x8632::misc-data-offset (% y)) (% imm0)) 45 (movl (@ x (% esp)) (% y)) 46 (movl ($ 0) (% temp0)) 47 @loop 48 (movl (@ x8632::misc-data-offset (% y) (% temp0)) (% temp1)) 49 (movl (% temp1) (@ (% imm0) (% temp0))) 50 (addl ($ 4) (% temp0)) 51 (cmpl (% temp0) (% l)) 52 (jnz @loop) 53 (mark-as-node temp1) 54 (single-value-return 3)) 55 56(CCL::defx8632lapfunction gmp-bignum-copy-negate-from-lisp 57 ((x 4) #|(ra 0)|# (y arg_y) (l arg_z)) 58 (mark-as-imm temp1) 59 (movl (@ x8632::misc-data-offset (% y)) (% imm0)) 60 (movl (@ x (% esp)) (% y)) 61 (movl ($ 0) (% temp0)) 62 @loop1 63 (movl (@ x8632::misc-data-offset (% y) (% temp0)) (% temp1)) 64 (cmpl ($ 0) (% temp1)) 65 (jnz @negate) 66 (movl ($ 0) (@ (% imm0) (% temp0))) 67 (addl ($ 4) (% temp0)) 68 (cmpl (% temp0) (% l)) 69 (jnz @loop1) 70 (jmp @return) 71 @negate 72 (neg (% temp1)) 73 (movl (% temp1) (@ (% imm0) (% temp0))) 74 (addl ($ 4) (% temp0)) 75 (cmpl (% temp0) (% l)) 76 (jz @return) 77 @loop2 78 (movl (@ x8632::misc-data-offset (% y) (% temp0)) (% temp1)) 79 (notl (% temp1)) 80 (movl (% temp1) (@ (% imm0) (% temp0))) 81 (addl ($ 4) (% temp0)) 82 (cmpl (% temp0) (% l)) 83 (jnz @loop2) 84 @return 85 (mark-as-node temp1) 86 (single-value-return 3)) 87 88(CCL::defx8632lapfunction gmp-bignum-copy-to-lisp 89 ((x 4) #|(ra 0)|# (y arg_y) (l arg_z)) 90 (mark-as-imm temp1) 91 (movl (@ x (% esp)) (% temp0)) 92 (movl (@ x8632::misc-data-offset (% temp0)) (% imm0)) 93 (movl ($ 0) (% temp0)) 94 @loop 95 (movl (@ (% imm0) (% temp0)) (% temp1)) 96 (movl (% temp1) (@ x8632::misc-data-offset (% y) (% temp0))) 97 (addl ($ 4) (% temp0)) 98 (cmpl (% temp0) (% l)) 99 (jnz @loop) 100 (mark-as-node temp1) 101 (single-value-return 3)) 102 103(CCL::defx8632lapfunction gmp-bignum-copy-negate-to-lisp 104 ((x 4) #|(ra 0)|# (y arg_y) (l arg_z)) 105 (mark-as-imm temp1) 106 (movl (@ x (% esp)) (% temp0)) 107 (movl (@ x8632::misc-data-offset (% temp0)) (% imm0)) 108 (movl ($ 0) (% temp0)) 109 @loop1 110 (movl (@ (% imm0) (% temp0)) (% temp1)) 111 (cmpl ($ 0) (% temp1)) 112 (jnz @negate) 113 (movl ($ 0) (@ x8632::misc-data-offset (% y) (% temp0))) 114 (addl ($ 4) (% temp0)) 115 (cmpl (% temp0) (% l)) 116 (jnz @loop1) 117 (jmp @return) 118 @negate 119 (neg (% temp1)) 120 (movl (% temp1) (@ x8632::misc-data-offset (% y) (% temp0))) 121 (addl ($ 4) (% temp0)) 122 (cmpl (% temp0) (% l)) 123 (jz @return) 124 @loop2 125 (movl (@ (% imm0) (% temp0)) (% temp1)) 126 (notl (% temp1)) 127 (movl (% temp1) (@ x8632::misc-data-offset (% y) (% temp0))) 128 (addl ($ 4) (% temp0)) 129 (cmpl (% temp0) (% l)) 130 (jnz @loop2) 131 @return 132 (mark-as-node temp1) 133 (single-value-return 3)) 134 135(CCL::defx8632lapfunction gmp-bignum-copy 136 ((x 4) #|(ra 0)|# (y arg_y) (l arg_z)) 137 (movl (@ x (% esp)) (% temp0)) 138 (movl ($ 0) (% temp1)) 139 @loop 140 (movl (@ x8632::misc-data-offset (% temp0) (% temp1)) (% imm0)) 141 (movl (% imm0) (@ x8632::misc-data-offset (% y) (% temp1))) 142 (addl ($ 4) (% temp1)) 143 (cmpl (% temp1) (% l)) 144 (jnz @loop) 145 (single-value-return 3)) 146 147) 148 149#+:X8664-TARGET 150(progn 151(CCL::defx86lapfunction gmp-bignum-copy-from-lisp 152 ((x arg_x) (y arg_y) (l arg_z)) 153 (movq ($ 0) (% imm1)) 154 (movq (@ x8664::misc-data-offset (% y)) (% imm2)) 155 @loop 156 (movq (@ x8664::misc-data-offset (% x) (% imm1)) (% imm0)) 157 (movq (% imm0) (@ (% imm2) (% imm1))) 158 (addq ($ 8) (% imm1)) 159 (cmpq (% imm1) (% l)) 160 (jnz @loop) 161 (single-value-return)) 162 163(CCL::defx86lapfunction gmp-bignum-copy-negate-from-lisp 164 ((x arg_x) (y arg_y) (l arg_z)) 165 (movq ($ 8) (% temp0)) 166 (andq (% l) (% temp0)) 167 (xorq (% temp0) (% l)) 168 (shrq ($ 1) (% l)) 169 (xorq (% imm1) (% imm1)) 170 (movq (@ x8664::misc-data-offset (% y)) (% imm2)) 171 (clc) 172 @loop1 173 (movq (@ x8664::misc-data-offset (% x) (% imm1)) (% imm0)) 174 (cmpq ($ 0) (% imm0)) 175 (jnz @negate) 176 (movq ($ 0) (@ (% imm2) (% imm1))) 177 (addq ($ 8) (% imm1)) 178 (cmpq (% imm1) (% l)) 179 (jnz @loop1) 180 (cmpq ($ 0) (% temp0)) 181 (jz @return) 182 (movslq (@ x8664::misc-data-offset (% x) (% imm1)) (% imm0)) 183 (xorq (% temp0) (% temp0)) 184 @negate 185 (neg (% imm0)) 186 (movq (% imm0) (@ (% imm2) (% imm1))) 187 (addq ($ 8) (% imm1)) 188 (cmpq (% imm1) (% l)) 189 (jle @finish) 190 @loop2 191 (movq (@ x8664::misc-data-offset (% x) (% imm1)) (% imm0)) 192 (notq (% imm0)) 193 (movq (% imm0) (@ (% imm2) (% imm1))) 194 (addq ($ 8) (% imm1)) 195 (cmpq (% imm1) (% l)) 196 (jnz @loop2) 197 @finish 198 (cmpq ($ 0) (% temp0)) 199 (jz @return) 200 (movslq (@ x8664::misc-data-offset (% x) (% imm1)) (% imm0)) 201 (notq (% imm0)) 202 (movq (% imm0) (@ (% imm2) (% imm1))) 203 @return 204 (single-value-return)) 205 206(CCL::defx86lapfunction gmp-bignum-copy-to-lisp 207 ((x arg_x) (y arg_y) (l arg_z)) 208 (movq ($ 0) (% imm1)) 209 (movq (@ x8664::misc-data-offset (% x)) (% imm2)) 210 @loop 211 (movq (@ (% imm2) (% imm1)) (% imm0)) 212 (movq (% imm0) (@ x8664::misc-data-offset (% y) (% imm1))) 213 (addq ($ 8) (% imm1)) 214 (cmpq (% imm1) (% l)) 215 (jnz @loop) 216 (single-value-return)) 217 218(CCL::defx86lapfunction gmp-bignum-copy-negate-to-lisp 219 ((x arg_x) (y arg_y) (l arg_z)) 220 (xorq (% imm1) (% imm1)) 221 (movq (@ x8664::misc-data-offset (% x)) (% imm2)) 222 @loop1 223 (movq (@ (% imm2) (% imm1)) (% imm0)) 224 (cmpq ($ 0) (% imm0)) 225 (jnz @negate) 226 (movq ($ 0) (@ x8664::misc-data-offset (% y) (% imm1))) 227 (addq ($ 8) (% imm1)) 228 (cmpq (% imm1) (% l)) 229 (jnz @loop1) 230 (jmp @return) 231 @negate 232 (neg (% imm0)) 233 (movq (% imm0) (@ x8664::misc-data-offset (% y) (% imm1))) 234 (addq ($ 8) (% imm1)) 235 (cmpq (% imm1) (% l)) 236 (jz @return) 237 @loop2 238 (movq (@ (% imm2) (% imm1)) (% imm0)) 239 (notq (% imm0)) 240 (movq (% imm0) (@ x8664::misc-data-offset (% y) (% imm1))) 241 (addq ($ 8) (% imm1)) 242 (cmpq (% imm1) (% l)) 243 (jnz @loop2) 244 @return 245 (single-value-return)) 246 247(CCL::defx86lapfunction gmp-bignum-copy 248 ((x arg_x) (y arg_y) (l arg_z)) 249 (movq ($ 0) (% imm1)) 250 @loop 251 (movq (@ x8664::misc-data-offset (% x) (% imm1)) (% imm0)) 252 (movq (% imm0) (@ x8664::misc-data-offset (% y) (% imm1))) 253 (addq ($ 8) (% imm1)) 254 (cmpq (% imm1) (% l)) 255 (jnz @loop) 256 (single-value-return)) 257) 258 259#+:ARM-TARGET 260(progn 261 262(CCL::defarmlapfunction gmp-bignum-copy-from-lisp 263 ((x arg_x) (y arg_y) (l arg_z)) 264 (mov imm1 (:$ 0)) 265 (ldr imm2 (:@ y (:$ arm::misc-data-offset))) 266 @loop 267 (add imm0 imm1 (:$ arm::misc-data-offset)) 268 (ldr imm0 (:@ x imm0)) 269 (str imm0 (:@ imm2 imm1)) 270 (add imm1 imm1 (:$ 4)) 271 (cmp imm1 l) 272 (bne @loop) 273 (bx lr)) 274 275(CCL::defarmlapfunction gmp-bignum-copy-negate-from-lisp 276 ((x arg_x) (y arg_y) (l arg_z)) 277 (mov imm1 (:$ 0)) 278 (ldr imm2 (:@ y (:$ arm::misc-data-offset))) 279 @loop1 280 (add imm0 imm1 (:$ arm::misc-data-offset)) 281 (ldr imm0 (:@ x imm0)) 282 (cmp imm0 (:$ 0)) 283 (bne @negate) 284 (str imm0 (:@ imm2 imm1)) 285 (add imm1 imm1 (:$ 4)) 286 (cmp imm1 l) 287 (bne @loop1) 288 (bx lr) 289 @negate 290 (rsb imm0 imm0 (:$ 0)) 291 (str imm0 (:@ imm2 imm1)) 292 (add imm1 imm1 (:$ 4)) 293 (cmp imm1 l) 294 (bne @loop2) 295 (bx lr) 296 @loop2 297 (add imm0 imm1 (:$ arm::misc-data-offset)) 298 (ldr imm0 (:@ x imm0)) 299 (mvn imm0 imm0) 300 (str imm0 (:@ imm2 imm1)) 301 (add imm1 imm1 (:$ 4)) 302 (cmp imm1 l) 303 (bne @loop2) 304 (bx lr)) 305 306(CCL::defarmlapfunction gmp-bignum-copy-to-lisp 307 ((x arg_x) (y arg_y) (l arg_z)) 308 (ldr imm2 (:@ x (:$ arm::misc-data-offset))) 309 @loop 310 (add l l (:$ -4)) 311 (ldr imm0 (:@ imm2 l)) 312 (add imm1 l (:$ arm::misc-data-offset)) 313 (str imm0 (:@ y imm1)) 314 (cmp l (:$ 0)) 315 (bne @loop) 316 (bx lr)) 317 318(CCL::defarmlapfunction gmp-bignum-copy-negate-to-lisp 319 ((x arg_x) (y arg_y) (l arg_z)) 320 (mov temp0 (:$ 0)) 321 (ldr imm2 (:@ x (:$ arm::misc-data-offset))) 322 @loop1 323 (ldr imm0 (:@ imm2 temp0)) 324 (cmp imm0 (:$ 0)) 325 (bne @negate) 326 (add imm1 temp0 (:$ arm::misc-data-offset)) 327 (str imm0 (:@ y imm1)) 328 (add temp0 temp0 (:$ 4)) 329 (cmp temp0 l) 330 (bne @loop1) 331 (bx lr) 332 @negate 333 (rsb imm0 imm0 (:$ 0)) 334 (add imm1 temp0 (:$ arm::misc-data-offset)) 335 (str imm0 (:@ y imm1)) 336 (add temp0 temp0 (:$ 4)) 337 (cmp temp0 l) 338 (bne @loop2) 339 (bx lr) 340 @loop2 341 (ldr imm0 (:@ imm2 temp0)) 342 (mvn imm0 imm0) 343 (add imm1 temp0 (:$ arm::misc-data-offset)) 344 (str imm0 (:@ y imm1)) 345 (add temp0 temp0 (:$ 4)) 346 (cmp temp0 l) 347 (bne @loop2) 348 (bx lr)) 349 350(CCL::defarmlapfunction gmp-bignum-copy 351 ((x arg_x) (y arg_y) (l arg_z)) 352 (mov imm1 (:$ arm::misc-data-offset)) 353 (add imm2 l (:$ arm::misc-data-offset)) 354 @loop 355 (ldr imm0 (:@ x imm1)) 356 (str imm0 (:@ y imm1)) 357 (add imm1 imm1 (:$ 4)) 358 (cmp imm1 imm2) 359 (bne @loop) 360 (bx lr)) 361) 362 363(defun gmp-bignum-isqrt (x) 364 (let* ((xl (digits_to_words (ccl::%bignum-length x))) 365 (rl (ceiling (+ 1 xl) 2)) 366 (xlb (words_to_bytes xl)) 367 (rlb (words_to_bytes rl)) 368 (rl2 (words_to_digits rl)) 369 (res (ccl::%alloc-misc rl2 (bignum_subtag)))) 370 (declare (type fixnum xl rl xlb rlb rl2)) 371 (ccl::%stack-block ((tx xlb) 372 (tr rlb)) 373 (gmp-bignum-copy-from-lisp x tx xl) 374 (ccl::external-call "gmp_wrap_isqrt" 375 :address tr :long rl 376 :address tx :long xl) 377 (gmp-bignum-copy-to-lisp tr res rl) 378 (ccl::%normalize-bignum-2 t res)))) 379 380(if (not (fboundp 'orig-multiply-bignums)) 381 (setf (symbol-function 'orig-multiply-bignums) 382 (symbol-function 'ccl::multiply-bignums))) 383 384(if (not (fboundp 'orig-bignum-gcd)) 385 (setf (symbol-function 'orig-positive-bignum-gcd) 386 (symbol-function 'ccl::%positive-bignum-bignum-gcd))) 387 388(if (not (fboundp 'orig-bignum-truncate)) 389 (setf (symbol-function 'orig-bignum-truncate) 390 (symbol-function 'ccl::bignum-truncate))) 391 392(if (not (fboundp 'orig-isqrt)) 393 (setf (symbol-function 'orig-isqrt) 394 (symbol-function 'common-lisp:isqrt))) 395 396(defun gmp-multiply-bignums(x y &optional res) 397 (let ((xl0 (ccl::%bignum-length x)) 398 (yl0 (ccl::%bignum-length y))) 399 (declare (type fixnum xl0 yl0)) 400 (if (< xl0 30) 401 (return-from gmp-multiply-bignums 402 (maybe_call_with_res_2 orig-multiply-bignums x y))) 403 (if (< yl0 30) 404 (return-from gmp-multiply-bignums 405 (maybe_call_with_res_2 orig-multiply-bignums y x))) 406 (if (< (+ xl0 yl0) 120) 407 (return-from gmp-multiply-bignums 408 (maybe_call_with_res_2 orig-multiply-bignums x y))) 409 (let* ((xl (digits_to_words xl0)) 410 (yl (digits_to_words yl0)) 411 (x-plusp nil) 412 (y-plusp nil) 413 (rl (+ xl yl)) 414 (rl2 (words_to_digits rl)) 415 (xlb 0) 416 (ylb 0) 417 (itmp 0) 418 (tmp nil) 419 (rlb 0) 420 (res (alloc_bignum rl2 res))) 421 (declare (type fixnum xl yl rl rl2 xlb ylb rlb itmp)) 422 ;;; XXX Does not work 423 ;;; (declare (dynamic-extent res)) 424 (if (< xl yl) 425 (progn 426 (setf itmp xl) 427 (setf xl yl) 428 (setf yl itmp) 429 (setf itmp xl0) 430 (setf xl0 yl0) 431 (setf yl0 itmp) 432 (setf tmp x) 433 (setf x y) 434 (setf y tmp))) 435 (setf xlb (words_to_bytes xl)) 436 (setf ylb (words_to_bytes yl)) 437 (setf rlb (+ xlb ylb)) 438 (setf x-plusp (is_plus x xl0)) 439 (setf y-plusp (is_plus y yl0)) 440 (ccl::%stack-block ((tx xlb) 441 (ty ylb) 442 (tr rlb)) 443 (if x-plusp 444 (gmp-bignum-copy-from-lisp x tx xl) 445 (gmp-bignum-copy-negate-from-lisp x tx xl0)) 446 (if y-plusp 447 (gmp-bignum-copy-from-lisp y ty yl) 448 (gmp-bignum-copy-negate-from-lisp y ty yl0)) 449 (ccl::external-call "__gmpn_mul" 450 :address tr 451 :address tx :long xl 452 :address ty :long yl) 453 (if (eq x-plusp y-plusp) 454 (gmp-bignum-copy-to-lisp tr res rl) 455 (gmp-bignum-copy-negate-to-lisp tr res rl)) 456 (ccl::%normalize-bignum-2 t res))))) 457 458 459(defun gmp-positive-bignum-gcd(x y &optional res) 460 (let* ((xl (digits_to_words (ccl::%bignum-length x))) 461 (yl (digits_to_words (ccl::%bignum-length y))) 462 (rl (if (< xl yl) xl yl)) 463 (xlb (words_to_bytes xl)) 464 (ylb (words_to_bytes yl)) 465 (rlb (+ xlb ylb)) 466 ) 467 (declare (type fixnum xl yl rl xlb ylb rlb)) 468 ;;; XXX Does not work 469 ;;; (declare (dynamic-extent res)) 470 (ccl::%stack-block ((tx xlb) 471 (ty ylb) 472 (tr rlb)) 473 (gmp-bignum-copy-from-lisp x tx xl) 474 (gmp-bignum-copy-from-lisp y ty yl) 475 (setf rl (ccl::external-call "gmp_wrap_gcd" 476 :address tr 477 :address tx :long xl 478 :address ty :long yl 479 :long)) 480 (setf res (alloc_bignum (words_to_digits rl) res)) 481 (gmp-bignum-copy-to-lisp tr res rl) 482 (ccl::%normalize-bignum-2 t res)))) 483;;; Tests 484;;; (truncate -51520943106947801344 17339521378867071488) 485;;; 486(defun gmp-bignum-truncate (x y &optional norem) 487 (declare (ignore norem)) 488 (if (and (eql (ccl::%bignum-length y) 1) 489 (eql (ccl::%typed-miscref :bignum y 0) 0)) 490 (error (make-condition 'division-by-zero 491 :operation 'gmp-bignum-truncate 492 :operands (list x 0)))) 493 (let* ((x-plusp (is_plus x (ccl::%bignum-length x))) 494 (y-plusp (is_plus y (ccl::%bignum-length y))) 495 (x (if x-plusp x (ccl::negate-bignum x nil))) 496 (y (if y-plusp y (ccl::negate-bignum y nil))) 497 (yl0 (ccl::%bignum-length y)) 498 (xl (digits_to_words (ccl::%bignum-length x))) 499 (yl2 (if (eq 0 (ccl::%typed-miscref :bignum y (- yl0 1))) 500 (- yl0 1) 501 yl0)) 502 (yl (digits_to_words yl2)) 503 (ql (+ 1 (- xl yl))) 504 (q nil) 505 (r nil) 506 (xlb (words_to_bytes xl)) 507 (ylb (words_to_bytes yl)) 508 (qlb (words_to_bytes ql))) 509 (declare (type fixnum xl yl yl0 yl2 ql xlb ylb qlb)) 510 (if (plusp (ccl::bignum-compare y x)) 511 (progn 512 (setf r (ccl::%alloc-misc (words_to_digits xl) (bignum_subtag))) 513 (gmp-bignum-copy x r xl) 514 (setf q 0)) 515 (ccl::%stack-block ((tx xlb) 516 (ty ylb) 517 (tq qlb) 518 (tr ylb)) 519 (gmp-bignum-copy-from-lisp x tx xl) 520 (gmp-bignum-copy-from-lisp y ty yl) 521 (ccl::external-call "__gmpn_tdiv_qr" 522 :address tq :address tr :long 0 523 :address tx :long xl 524 :address ty :long yl) 525 (setf q (ccl::%alloc-misc (words_to_digits ql) (bignum_subtag))) 526 (setf r (ccl::%alloc-misc yl0 (bignum_subtag))) 527 (gmp-bignum-copy-to-lisp tq q ql) 528 (gmp-bignum-copy-to-lisp tr r yl) 529 (if (> yl0 yl2) 530 (setf (ccl::%typed-miscref :bignum r (- yl0 1)) 0)) 531 (setf q (ccl::%normalize-bignum-2 t q)) 532 (setf r (ccl::%normalize-bignum-2 t r)))) 533 (let ((quotient (cond ((eq x-plusp y-plusp) q) 534 ((typep q 'fixnum) (the fixnum (- q))) 535 (t (ccl::negate-bignum-in-place q)))) 536 (rem (cond (x-plusp r) 537 ((typep r 'fixnum) (the fixnum (- r))) 538 (t (ccl::negate-bignum-in-place r))))) 539 (values (if (typep quotient 'fixnum) 540 quotient 541 (ccl::%normalize-bignum-2 t quotient)) 542 (if (typep rem 'fixnum) 543 rem 544 (ccl::%normalize-bignum-2 t rem)))))) 545 546) 547 548;;; Tests 549;;; (gmp-bignum-isqrt (expt 10 50)) 550;;; (gmp-bignum-isqrt (expt 2 127)) 551#+:sbcl 552(defun gmp-bignum-isqrt (x) 553 (let* ((len-x (sb-bignum::%bignum-length x)) 554 (len-res (ceiling (+ 1 len-x) 2)) 555 (res (sb-bignum::%allocate-bignum len-res))) 556 (declare (type fixnum len-x len-res)) 557 (sb-sys:with-pinned-objects (x res) 558 (let* ((addrx (sb-kernel:get-lisp-obj-address x)) 559 (sapx (sb-sys:int-sap addrx)) 560 (addr-res (sb-kernel:get-lisp-obj-address res)) 561 (sapr (sb-sys:int-sap addr-res))) 562 (sb-alien:alien-funcall 563 (sb-alien:extern-alien "gmp_wrap_sb_isqrt" 564 (sb-alien:function sb-alien:void (* t) (* t))) 565 sapx sapr))) 566 (sb-bignum::%normalize-bignum res len-res))) 567 568(defun gmp-isqrt (n) 569 (check-type n unsigned-byte) 570 (if (not #+:sbcl(sb-int:fixnump n) 571 #+:openmcl(ccl:fixnump n)) 572 (return-from gmp-isqrt (gmp-bignum-isqrt n))) 573 (locally 574 (declare (type fixnum n)) 575 (if (<= n 24) 576 (cond ((> n 15) 4) 577 ((> n 8) 3) 578 ((> n 3) 2) 579 ((> n 0) 1) 580 (t 0)) 581 (let* ((n-len-quarter (ash (integer-length n) -2)) 582 (n-half (ash n (- (ash n-len-quarter 1)))) 583 (n-half-isqrt (isqrt n-half)) 584 (init-value (ash (1+ n-half-isqrt) n-len-quarter))) 585 (declare (type fixnum n-len-quarter n-half 586 n-half-isqrt init-value)) 587 (loop 588 (let ((iterated-value 589 (ash (+ init-value (truncate n init-value)) -1))) 590 (unless (< iterated-value init-value) 591 (return init-value)) 592 (setq init-value iterated-value))))))) 593 594(defparameter *gmp-multiplication-initialized* nil) 595 596#+:openmcl 597(progn 598 599(defun install-gmp-multiplication() 600 (let ((package *PACKAGE*)) 601 (ccl:SET-DEVELOPMENT-ENVIRONMENT T) 602 (setf (symbol-function 'ccl::multiply-bignums) 603 (symbol-function 'gmp-multiply-bignums)) 604 (setf (symbol-function 'ccl::bignum-truncate) 605 (symbol-function 'gmp-bignum-truncate)) 606 (setf (symbol-function 'ccl::%positive-bignum-bignum-gcd) 607 (symbol-function 'gmp-positive-bignum-gcd)) 608 (setf (symbol-function 'common-lisp:isqrt) 609 (symbol-function 'gmp-isqrt)) 610 (ccl:SET-USER-ENVIRONMENT T) 611 (setf *PACKAGE* package)) 612) 613 614(defun uninstall-gmp-multiplication() 615 (let ((package *PACKAGE*)) 616 (ccl:SET-DEVELOPMENT-ENVIRONMENT T) 617 (setf (symbol-function 'ccl::multiply-bignums) 618 (symbol-function 'orig-multiply-bignums)) 619 (setf (symbol-function 'ccl::bignum-truncate) 620 (symbol-function 'orig-bignum-truncate)) 621 (setf (symbol-function 'ccl::%positive-bignum-bignum-gcd) 622 (symbol-function 'orig-positive-bignum-gcd)) 623 (setf (symbol-function 'common-lisp:isqrt) 624 (symbol-function 'orig-isqrt)) 625 (ccl:SET-USER-ENVIRONMENT T) 626 (setf *PACKAGE* package)) 627) 628 629) 630 631#+:sbcl 632(progn 633(if (not (fboundp 'orig-multiply-bignums)) 634 (setf (symbol-function 'orig-multiply-bignums) 635 (symbol-function 'sb-bignum::multiply-bignums))) 636 637(if (not (fboundp 'orig-bignum-gcd)) 638 (setf (symbol-function 'orig-bignum-gcd) 639 (symbol-function 'sb-bignum::bignum-gcd))) 640 641(if (not (fboundp 'orig-bignum-truncate)) 642 (setf (symbol-function 'orig-bignum-truncate) 643 (symbol-function 'sb-bignum::bignum-truncate))) 644 645(if (not (fboundp 'orig-isqrt)) 646 (setf (symbol-function 'orig-isqrt) 647 (symbol-function 'common-lisp:isqrt))) 648 649(defun gmp-multiply-bignums0 (a b) 650 ;;; (declare (type bignum-type a b)) 651 (let* ((a-plusp (sb-bignum::%bignum-0-or-plusp a 652 (sb-bignum::%bignum-length a))) 653 (b-plusp (sb-bignum::%bignum-0-or-plusp b 654 (sb-bignum::%bignum-length b))) 655 (a (if a-plusp a (sb-bignum::negate-bignum a))) 656 (b (if b-plusp b (sb-bignum::negate-bignum b))) 657 (len-a (sb-bignum::%bignum-length a)) 658 (len-b (sb-bignum::%bignum-length b)) 659 (len-res (+ len-a len-b)) 660 (res (sb-bignum::%allocate-bignum len-res)) 661 (negate-res (not (eql a-plusp b-plusp)))) 662 ;;; (declare (type bignum-index len-a len-b len-res)) 663 (if (< len-a len-b) 664 (let ((tmp a)) 665 (setf a b) 666 (setf b tmp))) 667 (sb-sys:with-pinned-objects (a b res) 668 (let* ((addra (sb-kernel:get-lisp-obj-address a)) 669 (sapa (sb-sys:int-sap addra)) 670 (addrb (sb-kernel:get-lisp-obj-address b)) 671 (sapb (sb-sys:int-sap addrb)) 672 (addr-res (sb-kernel:get-lisp-obj-address res)) 673 (sap-res (sb-sys:int-sap addr-res))) 674 (sb-alien:alien-funcall 675 (sb-alien:extern-alien "gmp_wrap_sb_mul" 676 (sb-alien:function sb-alien:void (* t) (* t) (* t))) 677 sapa sapb sap-res))) 678 (when negate-res (sb-bignum::negate-bignum-in-place res)) 679 (sb-bignum::%normalize-bignum res len-res) 680 681 ) 682) 683 684(defun gmp-multiply-bignums(x y) 685 (let ((xl0 (sb-bignum::%bignum-length x)) 686 (yl0 (sb-bignum::%bignum-length y))) 687 (declare (type fixnum xl0 yl0)) 688 (if (< xl0 10) 689 (return-from gmp-multiply-bignums (orig-multiply-bignums x y))) 690 (if (< yl0 10) 691 (return-from gmp-multiply-bignums (orig-multiply-bignums y x))) 692 (if (< (+ xl0 yl0) 40) 693 (return-from gmp-multiply-bignums (orig-multiply-bignums x y))) 694 (gmp-multiply-bignums0 x y))) 695 696(defun gmp-bignum-gcd(x y) 697 (let* ( 698 (nx (if (sb-bignum::%bignum-0-or-plusp x (sb-bignum::%bignum-length x)) 699 (sb-bignum::copy-bignum x) 700 (sb-bignum::negate-bignum x nil))) 701 (ny (if (sb-bignum::%bignum-0-or-plusp y (sb-bignum::%bignum-length y)) 702 (sb-bignum::copy-bignum y) 703 (sb-bignum::negate-bignum y nil))) 704 (xl (sb-bignum::%bignum-length nx)) 705 (yl (sb-bignum::%bignum-length ny)) 706 (rl (if (< xl yl) xl yl)) 707 (res (sb-bignum::%allocate-bignum rl))) 708 (sb-sys:with-pinned-objects (nx ny res) 709 (let* ((addrx (sb-kernel:get-lisp-obj-address nx)) 710 (sapx (sb-sys:int-sap addrx)) 711 (addry (sb-kernel:get-lisp-obj-address ny)) 712 (sapy (sb-sys:int-sap addry)) 713 (addr-res (sb-kernel:get-lisp-obj-address res)) 714 (sap-res (sb-sys:int-sap addr-res))) 715 (sb-alien:alien-funcall 716 (sb-alien:extern-alien "gmp_wrap_sb_gcd" 717 (sb-alien:function sb-alien:void (* t) (* t) (* t))) 718 sapx sapy sap-res))) 719 (sb-bignum::%normalize-bignum res rl) 720 ) 721) 722 723(defun test-bignum-gcd(x y) 724 (let ((res1 (orig-bignum-gcd x y)) 725 (res2 (gmp-bignum-gcd x y))) 726 (if (not (equal res1 res2)) 727 (format t 728 "Different results from gcd ~S ~S, orig ~S, gmp ~S ~%" 729 x y res1 res2)) 730 res2)) 731 732 733(defun gmp-bignum-truncate(x y) 734 (let* ( 735 (x-plusp (sb-bignum::%bignum-0-or-plusp x (sb-bignum::%bignum-length x))) 736 (y-plusp (sb-bignum::%bignum-0-or-plusp y (sb-bignum::%bignum-length y))) 737 (nx (if x-plusp x 738 (sb-bignum::negate-bignum x nil))) 739 (ny (if y-plusp y 740 (sb-bignum::negate-bignum y nil))) 741 (len-x (sb-bignum::%bignum-length nx)) 742 (len-y (sb-bignum::%bignum-length ny)) 743 (q nil) 744 (r nil) 745 ) 746 (if (plusp (sb-bignum::bignum-compare ny nx)) 747 (progn 748 (setf q 0) 749 (setf r (if y-plusp (sb-bignum::copy-bignum nx) nx)) 750 ) 751 (let* ( 752 (len-r len-y) 753 (len-y (if (eql 0 (sb-bignum::%bignum-ref ny (- len-y 1))) 754 (- len-y 1) 755 len-y)) 756 (len-q (+ 1 (- len-x len-y))) 757 (nq (sb-bignum::%allocate-bignum len-q)) 758 (nr (sb-bignum::%allocate-bignum len-r))) 759 (sb-sys:with-pinned-objects (nx ny nq nr) 760 (let* ( 761 (addrx (sb-kernel:get-lisp-obj-address nx)) 762 (sapx (sb-sys:int-sap addrx)) 763 (addry (sb-kernel:get-lisp-obj-address ny)) 764 (sapy (sb-sys:int-sap addry)) 765 (addr-quo (sb-kernel:get-lisp-obj-address nq)) 766 (sapq (sb-sys:int-sap addr-quo)) 767 (addr-rem (sb-kernel:get-lisp-obj-address nr)) 768 (sapr (sb-sys:int-sap addr-rem))) 769 (sb-alien:alien-funcall 770 (sb-alien:extern-alien "gmp_wrap_sb_div_rem" 771 (sb-alien:function sb-alien:void (* t) (* t) (* t) (* t))) 772 sapx sapy sapq sapr))) 773 (setf q (sb-bignum::%normalize-bignum nq len-q)) 774 (setf r (sb-bignum::%normalize-bignum nr len-r)))) 775 (let ((quotient (cond ((eql x-plusp y-plusp) q) 776 ((typep q 'fixnum) (the fixnum (- q))) 777 (t (sb-bignum::negate-bignum-in-place q)))) 778 (rem (cond (x-plusp r) 779 ((typep r 'fixnum) (the fixnum (- r))) 780 (t (sb-bignum::negate-bignum-in-place r))))) 781 (values (if (typep quotient 'fixnum) 782 quotient 783 (sb-bignum::%normalize-bignum quotient 784 (sb-bignum::%bignum-length quotient))) 785 (if (typep rem 'fixnum) 786 rem 787 (sb-bignum::%normalize-bignum rem 788 (sb-bignum::%bignum-length rem))))))) 789 790 791 792#| 793;;; Tests 794;;; (truncate -51520943106947801344 17339521378867071488) 795;;; (truncate 23215968175662844254 12149601698671348626) 796;;; (truncate 1666974137583209287393566720 -2023369608) 797|# 798 799(defun install-gmp-multiplication() 800 (sb-ext:unlock-package "SB-BIGNUM") 801 (setf (symbol-function 'sb-bignum::multiply-bignums) 802 (symbol-function 'gmp-multiply-bignums)) 803 (setf (symbol-function 'sb-bignum::bignum-truncate) 804 (symbol-function 'gmp-bignum-truncate)) 805 (setf (symbol-function 'sb-bignum::bignum-gcd) 806 (symbol-function 'gmp-bignum-gcd)) 807 (sb-ext:lock-package "SB-BIGNUM") 808 (sb-ext:unlock-package "COMMON-LISP") 809 (setf (symbol-function 'common-lisp:isqrt) 810 (symbol-function 'gmp-isqrt)) 811 (sb-ext:lock-package "COMMON-LISP")) 812 813(defun uninstall-gmp-multiplication() 814 (sb-ext:unlock-package "SB-BIGNUM") 815 (setf (symbol-function 'sb-bignum::multiply-bignums) 816 (symbol-function 'orig-multiply-bignums)) 817 (setf (symbol-function 'sb-bignum::bignum-truncate) 818 (symbol-function 'orig-bignum-truncate)) 819 (setf (symbol-function 'sb-bignum::bignum-gcd) 820 (symbol-function 'orig-bignum-gcd)) 821 (sb-ext:lock-package "SB-BIGNUM") 822 (sb-ext:unlock-package "COMMON-LISP") 823 (setf (symbol-function 'common-lisp:isqrt) 824 (symbol-function 'orig-isqrt)) 825 (sb-ext:lock-package "COMMON-LISP"))) 826 827(defun init-gmp(wrapper-lib) 828 (if (not *gmp-multiplication-initialized*) 829 (if (ignore-errors (|quiet_load_alien| "libgmp.so") t) 830 (if (ignore-errors 831 (|quiet_load_alien| wrapper-lib) t) 832 (progn 833 (install-gmp-multiplication) 834 (setf *gmp-multiplication-initialized* t)))))) 835) 836