contrib/sb-mpfr/mpfr.lisp

   1 (defpackage :sb-mpfr
   2   (:use "COMMON-LISP" "SB-ALIEN" "SB-C-CALL")
   3   (:import-from "SB-GMP"
   4                 #:make-gmp-rstate
   5                 #:make-gmp-rstate-lc
   6                 #:rand-seed)
   7   (:export
   8    ;; +mpfr-precision+ is a pseudo-constant
   9    ;; Do not set via LET but use the two methods below!
  10    #:+mpfr-precision+
  11    #:set-precision
  12    #:with-precision
  13    ;; parameters and types
  14    #:*mpfr-rnd*
  15    #:*mpfr-base*
  16    #:mpfr-float
  17    ;; arithmetic operations
  18    #:make-mpfr-float
  19    #:mpfr-float-to-string
  20    #:add
  21    #:sub
  22    #:mul
  23    #:square
  24    #:div
  25    #:sqrt
  26    #:reciprocal-sqrt
  27    #:cubic-root
  28    #:k-root
  29    #:power
  30    #:negate
  31    #:abs
  32    #:dim
  33    #:mul-2-raised
  34    #:div-2-raised
  35    ;; special functions
  36    #:log
  37    #:log2
  38    #:log10
  39    #:exp
  40    #:exp2
  41    #:exp10
  42    #:cos
  43    #:sin
  44    #:tan
  45    #:sin-cos
  46    #:sec
  47    #:csc
  48    #:cot
  49    #:acos
  50    #:asin
  51    #:atan
  52    #:cosh
  53    #:sinh
  54    #:tanh
  55    #:sinh-cosh
  56    #:sech
  57    #:csch
  58    #:coth
  59    #:acosh
  60    #:asinh
  61    #:atanh
  62    #:fac
  63    #:log1p
  64    #:expm1
  65    #:eint
  66    #:li2
  67    #:gamma
  68    #:lngamma
  69    #:digamma
  70    #:zeta
  71    #:erf
  72    #:erfc
  73    #:j0
  74    #:j1
  75    #:jn
  76    #:y0
  77    #:y1
  78    #:yn
  79    #:ai
  80    #:arithmetic-geometric-mean
  81    #:hypot
  82    #:fma
  83    #:fms
  84    #:sum
  85    ;; comparison functions and predicates
  86    #:nan-p
  87    #:infinityp
  88    #:integerp
  89    #:numberp
  90    #:zerop
  91    #:regularp
  92    #:compare
  93    #:compare-2exp
  94    #:compare-abs
  95    #:>
  96    #:>=
  97    #:<
  98    #:<=
  99    #:=
 100    #:/=
 101    #:unorderedp
 102    ;; constants
 103    #:const-log2
 104    #:const-pi
 105    #:const-euler
 106    #:const-catalan
 107    ;; miscellaneous functions
 108    #:clear-underflow
 109    #:clear-overflow
 110    #:clear-div-by-zero
 111    #:clear-nan-flag
 112    #:clear-inex-flag
 113    #:clear-erange-flag
 114    #:set-underflow-flag
 115    #:set-overflow-flag
 116    #:set-div-by-zero-flag
 117    #:set-nan-flag
 118    #:set-inex-flag
 119    #:set-erange-flag
 120    #:clear-flags
 121    #:underflowp
 122    #:overflowp
 123    #:div-by-zero-p
 124    #:nanflag-p
 125    #:inexflag-p
 126    #:erangeflag-p
 127    ;; random number generation
 128    #:urandomb
 129    #:urandom
 130    ;; rounding
 131    #:rounded-int
 132    #:rounded-int-ceiling
 133    #:rounded-int-floor
 134    #:rounded-int-round
 135    #:rounded-int-truncate
 136    #:fractional
 137    #:ceiling
 138    #:floor
 139    #:round
 140    #:truncate
 141    #:modf
 142    #:fmod
 143    #:remainder
 144    #:remainder-quot
 145    ;; conversion
 146    #:coerce
 147    ;; special constants
 148    #:*mpfr-version*
 149    #:*mpfr-features*)
 150   (:shadow
 151    #:sqrt
 152    #:abs
 153    #:log
 154    #:exp
 155    #:cos
 156    #:sin
 157    #:tan
 158    #:acos
 159    #:asin
 160    #:atan
 161    #:cosh
 162    #:sinh
 163    #:tanh
 164    #:acosh
 165    #:asinh
 166    #:atanh
 167    #:numberp
 168    #:integerp
 169    #:zerop
 170    #:>
 171    #:>=
 172    #:<
 173    #:<=
 174    #:=
 175    #:/=
 176    #:ceiling
 177    #:floor
 178    #:round
 179    #:truncate
 180    #:coerce))
 181
 182 (in-package :sb-mpfr)
 183
 184 (defvar +mpfr-precision+)
 185 (defvar *mpfr-version* nil)
 186 (defvar *mpfr-features* nil)
 187
 188 (defun try-load-shared-object (pathname)
 189   (handler-case
 190       (load-shared-object pathname :dont-save t)
 191     (error (e)
 192       nil)))
 193
 194 (defun %load-mpfr ()
 195   (or (some #'try-load-shared-object
 196             #-(or win32 darwin) '("libmpfr.so" "libmpfr.so.4")
 197             #+darwin '("libmpfr.dylib" "libmpfr.4.dylib")
 198             #+win32 '("mpfr.dll"))
 199       (warn "MPFR not loaded.")))
 200
 201 (defun load-mpfr (&key (persistently t))
 202   (setf *mpfr-version* nil
 203         *mpfr-features* nil
 204         *features* (set-difference *features* '(:sb-mpfr)))
 205   (when persistently
 206     (pushnew 'load-mpfr sb-ext:*init-hooks*)
 207     ;(pushnew 'uninstall-mpfr-funs sb-ext:*save-hooks*)
 208     )
 209   (let ((success (%load-mpfr)))
 210     (when success
 211       (setf *mpfr-version*
 212             (alien-funcall (extern-alien "mpfr_get_version"
 213                                          (function c-string)))))
 214     (cond ((null *mpfr-version*))
 215           ((string<= *mpfr-version* "3.1")
 216            (warn "SB-MPFR requires at least MPFR version 3.1")
 217            (setf success nil))
 218           (t
 219            (setf +mpfr-precision+
 220                  (alien-funcall (extern-alien "mpfr_get_default_prec"
 221                                               (function long))))
 222            (pushnew :sb-mpfr *mpfr-features*)
 223            (setf *features* (union *features* *mpfr-features*))))
 224     success))
 225
 226 (load-mpfr)
 227
 228 ;;; types and initialization
 229
 230 (define-alien-type nil
 231     (struct mpfrfloat
 232             (mpfr_prec long)
 233             (mpfr_sign int)
 234             (mpfr_exp long)
 235             (mpfr_d (* unsigned-long))))
 236
 237 (define-alien-type mpfr_rnd_enum
 238   (enum mpfr_rnd
 239         (:mpfr_rndna -1)
 240         (:mpfr_rndn 0)
 241         (:mpfr_rndz 1)
 242         (:mpfr_rndu 2)
 243         (:mpfr_rndd 3)
 244         (:mpfr_rnda 4)
 245         (:mpfr_rndf 5)))
 246
 247 (declaim (inline mpfr_init2
 248                  mpfr_clear
 249                  mpfr_set
 250                  mpfr_set_ui
 251                  mpfr_set_si
 252                  mpfr_set_flt
 253                  mpfr_set_d
 254                  mpfr_set_z
 255                  mpfr_set_q
 256                  mpfr_set_nan
 257                  mpfr_set_inf
 258                  mpfr_set_zero
 259                  mpfr_set_default_prec
 260                  mpfr_get_str))
 261
 262 (define-alien-routine mpfr_init2 void
 263   (x (* (struct mpfrfloat)))
 264   (precision long))
 265
 266 (define-alien-routine mpfr_clear void
 267   (x (* (struct mpfrfloat))))
 268
 269 ;;; conversion functions
 270
 271 (define-alien-routine mpfr_set void
 272   (x (* (struct mpfrfloat)))
 273   (op (* (struct mpfrfloat)))
 274   (rnd mpfr_rnd_enum))
 275
 276 (define-alien-routine mpfr_set_ui void
 277   (x (* (struct mpfrfloat)))
 278   (op unsigned-long)
 279   (rnd mpfr_rnd_enum))
 280
 281 (define-alien-routine mpfr_set_si void
 282   (x (* (struct mpfrfloat)))
 283   (op long)
 284   (rnd mpfr_rnd_enum))
 285
 286 (define-alien-routine mpfr_set_flt void
 287   (x (* (struct mpfrfloat)))
 288   (op float)
 289   (rnd mpfr_rnd_enum))
 290
 291 (define-alien-routine mpfr_set_d void
 292   (x (* (struct mpfrfloat)))
 293   (op double-float)
 294   (rnd mpfr_rnd_enum))
 295
 296 (define-alien-routine mpfr_set_z void
 297   (x (* (struct mpfrfloat)))
 298   (op (* (struct sb-gmp::gmpint)))
 299   (rnd mpfr_rnd_enum))
 300
 301 (define-alien-routine mpfr_set_q void
 302   (x (* (struct mpfrfloat)))
 303   (op (* (struct sb-gmp::gmprat)))
 304   (rnd mpfr_rnd_enum))
 305
 306 (define-alien-routine mpfr_set_nan void
 307   (x (* (struct mpfrfloat))))
 308
 309 (define-alien-routine mpfr_set_inf void
 310   (x (* (struct mpfrfloat)))
 311   (sign int))
 312
 313 (define-alien-routine mpfr_set_zero void
 314   (x (* (struct mpfrfloat)))
 315   (sign int))
 316
 317 (define-alien-routine mpfr_set_default_prec void
 318   (prec long))
 319
 320 (define-alien-routine mpfr_get_str (* char)
 321   (str (* char))
 322   (exp (* long))
 323   (base int)
 324   (n unsigned)
 325   (x (* (struct mpfrfloat)))
 326   (rnd mpfr_rnd_enum))
 327
 328 ;;; conversion functions
 329
 330 (declaim (inline mpfr_get_flt
 331                  mpfr_get_d
 332                  mpfr_get_si
 333                  mpfr_get_ui
 334                  mpfr_get_z
 335                  mpfr_free_str
 336                  mpfr_set_str))
 337
 338 (define-alien-routine mpfr_get_flt float
 339   (op (* (struct mpfrfloat)))
 340   (rnd mpfr_rnd_enum))
 341
 342 (define-alien-routine mpfr_get_d double
 343   (op (* (struct mpfrfloat)))
 344   (rnd mpfr_rnd_enum))
 345
 346 (define-alien-routine mpfr_get_si long
 347   (op (* (struct mpfrfloat)))
 348   (rnd mpfr_rnd_enum))
 349
 350 (define-alien-routine mpfr_get_ui unsigned-long
 351   (op (* (struct mpfrfloat)))
 352   (rnd mpfr_rnd_enum))
 353
 354 (define-alien-routine mpfr_get_z int
 355   (res (* (struct sb-gmp::gmpint)))
 356   (op (* (struct mpfrfloat)))
 357   (rnd mpfr_rnd_enum))
 358
 359 (define-alien-routine mpfr_free_str void
 360   (str (* char)))
 361
 362 (define-alien-routine mpfr_set_str int
 363   (x (* (struct mpfrfloat)))
 364   (str c-string)
 365   (base int)
 366   (rnd mpfr_rnd_enum))
 367
 368 ;;; arithmetic functions
 369
 370 (declaim (inline mpfr_add
 371                  mpfr_add_ui
 372                  mpfr_add_si
 373                  mpfr_add_d
 374                  mpfr_add_z
 375                  mpfr_add_q
 376                  mpfr_sub
 377                  mpfr_ui_sub
 378                  mpfr_sub_ui
 379                  mpfr_si_sub
 380                  mpfr_sub_si
 381                  mpfr_d_sub
 382                  mpfr_sub_d
 383                  mpfr_z_sub
 384                  mpfr_sub_z
 385                  mpfr_sub_q
 386                  mpfr_mul
 387                  mpfr_mul_ui
 388                  mpfr_mul_si
 389                  mpfr_mul_d
 390                  mpfr_mul_z
 391                  mpfr_mul_q
 392                  mpfr_sqr
 393                  mpfr_div
 394                  mpfr_ui_div
 395                  mpfr_div_ui
 396                  mpfr_si_div
 397                  mpfr_div_si
 398                  mpfr_d_div
 399                  mpfr_div_d
 400                  mpfr_div_z
 401                  mpfr_div_q
 402                  mpfr_sqrt
 403                  mpfr_sqrt_ui
 404                  mpfr_rec_sqrt
 405                  mpfr_cbrt
 406                  mpfr_root
 407                  mpfr_pow
 408                  mpfr_pow_ui
 409                  mpfr_pow_si
 410                  mpfr_pow_z
 411                  mpfr_ui_pow
 412                  mpfr_neg
 413                  mpfr_abs
 414                  mpfr_dim
 415                  mpfr_mul_2ui
 416                  mpfr_mul_2si
 417                  mpfr_div_2ui
 418                  mpfr_div_2si))
 419
 420 (define-alien-routine mpfr_add int
 421   (r (* (struct mpfrfloat)))
 422   (op1 (* (struct mpfrfloat)))
 423   (op2 (* (struct mpfrfloat)))
 424   (rnd mpfr_rnd_enum))
 425
 426 (define-alien-routine mpfr_add_ui int
 427   (r (* (struct mpfrfloat)))
 428   (op1 (* (struct mpfrfloat)))
 429   (op2 unsigned-long)
 430   (rnd mpfr_rnd_enum))
 431
 432 (define-alien-routine mpfr_add_si int
 433   (r (* (struct mpfrfloat)))
 434   (op1 (* (struct mpfrfloat)))
 435   (op2 long)
 436   (rnd mpfr_rnd_enum))
 437
 438 (define-alien-routine mpfr_add_d int
 439   (r (* (struct mpfrfloat)))
 440   (op1 (* (struct mpfrfloat)))
 441   (op2 double-float)
 442   (rnd mpfr_rnd_enum))
 443
 444 (define-alien-routine mpfr_add_z int
 445   (r (* (struct mpfrfloat)))
 446   (op1 (* (struct mpfrfloat)))
 447   (op2 (* (struct sb-gmp::gmpint)))
 448   (rnd mpfr_rnd_enum))
 449
 450 (define-alien-routine mpfr_add_q int
 451   (r (* (struct mpfrfloat)))
 452   (op1 (* (struct mpfrfloat)))
 453   (op2 (* (struct sb-gmp::gmprat)))
 454   (rnd mpfr_rnd_enum))
 455
 456 (define-alien-routine mpfr_sub int
 457   (r (* (struct mpfrfloat)))
 458   (op1 (* (struct mpfrfloat)))
 459   (op2 (* (struct mpfrfloat)))
 460   (rnd mpfr_rnd_enum))
 461
 462 (define-alien-routine mpfr_ui_sub int
 463   (r (* (struct mpfrfloat)))
 464   (op1 unsigned-long)
 465   (op2 (* (struct mpfrfloat)))
 466   (rnd mpfr_rnd_enum))
 467
 468 (define-alien-routine mpfr_sub_ui int
 469   (r (* (struct mpfrfloat)))
 470   (op1 (* (struct mpfrfloat)))
 471   (op2 unsigned-long)
 472   (rnd mpfr_rnd_enum))
 473
 474 (define-alien-routine mpfr_si_sub int
 475   (r (* (struct mpfrfloat)))
 476   (op1 long)
 477   (op2 (* (struct mpfrfloat)))
 478   (rnd mpfr_rnd_enum))
 479
 480 (define-alien-routine mpfr_sub_si int
 481   (r (* (struct mpfrfloat)))
 482   (op1 (* (struct mpfrfloat)))
 483   (op2 long)
 484   (rnd mpfr_rnd_enum))
 485
 486 (define-alien-routine mpfr_d_sub int
 487   (r (* (struct mpfrfloat)))
 488   (op1 double-float)
 489   (op2 (* (struct mpfrfloat)))
 490   (rnd mpfr_rnd_enum))
 491
 492 (define-alien-routine mpfr_sub_d int
 493   (r (* (struct mpfrfloat)))
 494   (op1 (* (struct mpfrfloat)))
 495   (op2 double-float)
 496   (rnd mpfr_rnd_enum))
 497
 498 (define-alien-routine mpfr_z_sub int
 499   (r (* (struct mpfrfloat)))
 500   (op1 (* (struct sb-gmp::gmpint)))
 501   (op2 (* (struct mpfrfloat)))
 502   (rnd mpfr_rnd_enum))
 503
 504 (define-alien-routine mpfr_sub_z int
 505   (r (* (struct mpfrfloat)))
 506   (op1 (* (struct mpfrfloat)))
 507   (op2 (* (struct sb-gmp::gmpint)))
 508   (rnd mpfr_rnd_enum))
 509
 510 (define-alien-routine mpfr_sub_q int
 511   (r (* (struct mpfrfloat)))
 512   (op1 (* (struct mpfrfloat)))
 513   (op2 (* (struct sb-gmp::gmprat)))
 514   (rnd mpfr_rnd_enum))
 515
 516 (define-alien-routine mpfr_mul int
 517   (r (* (struct mpfrfloat)))
 518   (op1 (* (struct mpfrfloat)))
 519   (op2 (* (struct mpfrfloat)))
 520   (rnd mpfr_rnd_enum))
 521
 522 (define-alien-routine mpfr_mul_ui int
 523   (r (* (struct mpfrfloat)))
 524   (op1 (* (struct mpfrfloat)))
 525   (op2 unsigned-long)
 526   (rnd mpfr_rnd_enum))
 527
 528 (define-alien-routine mpfr_mul_si int
 529   (r (* (struct mpfrfloat)))
 530   (op1 (* (struct mpfrfloat)))
 531   (op2 long)
 532   (rnd mpfr_rnd_enum))
 533
 534 (define-alien-routine mpfr_mul_d int
 535   (r (* (struct mpfrfloat)))
 536   (op1 (* (struct mpfrfloat)))
 537   (op2 double-float)
 538   (rnd mpfr_rnd_enum))
 539
 540 (define-alien-routine mpfr_mul_z int
 541   (r (* (struct mpfrfloat)))
 542   (op1 (* (struct mpfrfloat)))
 543   (op2 (* (struct sb-gmp::gmpint)))
 544   (rnd mpfr_rnd_enum))
 545
 546 (define-alien-routine mpfr_mul_q int
 547   (r (* (struct mpfrfloat)))
 548   (op1 (* (struct mpfrfloat)))
 549   (op2 (* (struct sb-gmp::gmprat)))
 550   (rnd mpfr_rnd_enum))
 551
 552 (define-alien-routine mpfr_sqr int
 553   (r (* (struct mpfrfloat)))
 554   (op (* (struct mpfrfloat)))
 555   (rnd mpfr_rnd_enum))
 556
 557 (define-alien-routine mpfr_div int
 558   (r (* (struct mpfrfloat)))
 559   (op1 (* (struct mpfrfloat)))
 560   (op2 (* (struct mpfrfloat)))
 561   (rnd mpfr_rnd_enum))
 562
 563 (define-alien-routine mpfr_ui_div int
 564   (r (* (struct mpfrfloat)))
 565   (op1 unsigned-long)
 566   (op2 (* (struct mpfrfloat)))
 567   (rnd mpfr_rnd_enum))
 568
 569 (define-alien-routine mpfr_div_ui int
 570   (r (* (struct mpfrfloat)))
 571   (op1 (* (struct mpfrfloat)))
 572   (op2 unsigned-long)
 573   (rnd mpfr_rnd_enum))
 574
 575 (define-alien-routine mpfr_si_div int
 576   (r (* (struct mpfrfloat)))
 577   (op1 long)
 578   (op2 (* (struct mpfrfloat)))
 579   (rnd mpfr_rnd_enum))
 580
 581 (define-alien-routine mpfr_div_si int
 582   (r (* (struct mpfrfloat)))
 583   (op1 (* (struct mpfrfloat)))
 584   (op2 long)
 585   (rnd mpfr_rnd_enum))
 586
 587 (define-alien-routine mpfr_d_div int
 588   (r (* (struct mpfrfloat)))
 589   (op1 double-float)
 590   (op2 (* (struct mpfrfloat)))
 591   (rnd mpfr_rnd_enum))
 592
 593 (define-alien-routine mpfr_div_d int
 594   (r (* (struct mpfrfloat)))
 595   (op1 (* (struct mpfrfloat)))
 596   (op2 double-float)
 597   (rnd mpfr_rnd_enum))
 598
 599 (define-alien-routine mpfr_div_z int
 600   (r (* (struct mpfrfloat)))
 601   (op1 (* (struct mpfrfloat)))
 602   (op2 (* (struct sb-gmp::gmpint)))
 603   (rnd mpfr_rnd_enum))
 604
 605 (define-alien-routine mpfr_div_q int
 606   (r (* (struct mpfrfloat)))
 607   (op1 (* (struct mpfrfloat)))
 608   (op2 (* (struct sb-gmp::gmprat)))
 609   (rnd mpfr_rnd_enum))
 610
 611 (define-alien-routine mpfr_sqrt int
 612   (r (* (struct mpfrfloat)))
 613   (op (* (struct mpfrfloat)))
 614   (rnd mpfr_rnd_enum))
 615
 616 (define-alien-routine mpfr_sqrt_ui int
 617   (r (* (struct mpfrfloat)))
 618   (op unsigned-long)
 619   (rnd mpfr_rnd_enum))
 620
 621 (define-alien-routine mpfr_rec_sqrt int
 622   (r (* (struct mpfrfloat)))
 623   (op (* (struct mpfrfloat)))
 624   (rnd mpfr_rnd_enum))
 625
 626 (define-alien-routine mpfr_cbrt int
 627   (r (* (struct mpfrfloat)))
 628   (op (* (struct mpfrfloat)))
 629   (rnd mpfr_rnd_enum))
 630
 631 (define-alien-routine mpfr_root int
 632   (r (* (struct mpfrfloat)))
 633   (op (* (struct mpfrfloat)))
 634   (k unsigned-long)
 635   (rnd mpfr_rnd_enum))
 636
 637 (define-alien-routine mpfr_pow int
 638   (r (* (struct mpfrfloat)))
 639   (op1 (* (struct mpfrfloat)))
 640   (op2 (* (struct mpfrfloat)))
 641   (rnd mpfr_rnd_enum))
 642
 643 (define-alien-routine mpfr_pow_ui int
 644   (r (* (struct mpfrfloat)))
 645   (op1 (* (struct mpfrfloat)))
 646   (op2 unsigned-long)
 647   (rnd mpfr_rnd_enum))
 648
 649 (define-alien-routine mpfr_pow_si int
 650   (r (* (struct mpfrfloat)))
 651   (op1 (* (struct mpfrfloat)))
 652   (op2 long)
 653   (rnd mpfr_rnd_enum))
 654
 655 (define-alien-routine mpfr_pow_z int
 656   (r (* (struct mpfrfloat)))
 657   (op1 (* (struct mpfrfloat)))
 658   (op2 (* (struct sb-gmp::gmpint)))
 659   (rnd mpfr_rnd_enum))
 660
 661 (define-alien-routine mpfr_ui_pow int
 662   (r (* (struct mpfrfloat)))
 663   (op1 unsigned-long)
 664   (op2 (* (struct mpfrfloat)))
 665   (rnd mpfr_rnd_enum))
 666
 667 (define-alien-routine mpfr_neg int
 668   (r (* (struct mpfrfloat)))
 669   (op (* (struct mpfrfloat)))
 670   (rnd mpfr_rnd_enum))
 671
 672 (define-alien-routine mpfr_abs int
 673   (r (* (struct mpfrfloat)))
 674   (op (* (struct mpfrfloat)))
 675   (rnd mpfr_rnd_enum))
 676
 677 (define-alien-routine mpfr_dim int
 678   (r (* (struct mpfrfloat)))
 679   (op1 (* (struct mpfrfloat)))
 680   (op2 (* (struct mpfrfloat)))
 681   (rnd mpfr_rnd_enum))
 682
 683 (define-alien-routine mpfr_mul_2ui int
 684   (r (* (struct mpfrfloat)))
 685   (op1 (* (struct mpfrfloat)))
 686   (op2 unsigned-long)
 687   (rnd mpfr_rnd_enum))
 688
 689 (define-alien-routine mpfr_mul_2si int
 690   (r (* (struct mpfrfloat)))
 691   (op1 (* (struct mpfrfloat)))
 692   (op2 long)
 693   (rnd mpfr_rnd_enum))
 694
 695 (define-alien-routine mpfr_div_2ui int
 696   (r (* (struct mpfrfloat)))
 697   (op1 (* (struct mpfrfloat)))
 698   (op2 unsigned-long)
 699   (rnd mpfr_rnd_enum))
 700
 701 (define-alien-routine mpfr_div_2si int
 702   (r (* (struct mpfrfloat)))
 703   (op1 (* (struct mpfrfloat)))
 704   (op2 long)
 705   (rnd mpfr_rnd_enum))
 706
 707 ;;; special functions
 708
 709 (defmacro define-onearg-mpfr-int (funs)
 710   (loop for i in funs collect `(define-alien-routine ,i int
 711                                  (r (* (struct mpfrfloat)))
 712                                  (op (* (struct mpfrfloat)))
 713                                  (rnd mpfr_rnd_enum))
 714           into defines
 715         finally (return `(progn
 716                            (declaim (inline ,@funs))
 717                            ,@defines))))
 718
 719 (define-onearg-mpfr-int
 720     (mpfr_log
 721      mpfr_log2
 722      mpfr_log10
 723      mpfr_exp
 724      mpfr_exp2
 725      mpfr_exp10
 726      mpfr_cos
 727      mpfr_sin
 728      mpfr_tan
 729      mpfr_sec
 730      mpfr_csc
 731      mpfr_cot
 732      mpfr_acos
 733      mpfr_asin
 734      mpfr_atan
 735      mpfr_cosh
 736      mpfr_sinh
 737      mpfr_tanh
 738      mpfr_sech
 739      mpfr_csch
 740      mpfr_coth
 741      mpfr_acosh
 742      mpfr_asinh
 743      mpfr_atanh
 744      mpfr_log1p
 745      mpfr_expm1
 746      mpfr_eint
 747      mpfr_li2
 748      mpfr_gamma
 749      mpfr_lngamma
 750      mpfr_digamma
 751      mpfr_zeta
 752      mpfr_erf
 753      mpfr_erfc
 754      mpfr_j0
 755      mpfr_j1
 756      mpfr_y0
 757      mpfr_y1
 758      mpfr_ai))
 759
 760 (defmacro define-twoarg-mpfr-int (funs)
 761   (loop for i in funs collect `(define-alien-routine ,i int
 762                                  (r (* (struct mpfrfloat)))
 763                                  (op1 (* (struct mpfrfloat)))
 764                                  (op2 (* (struct mpfrfloat)))
 765                                  (rnd mpfr_rnd_enum))
 766           into defines
 767         finally (return `(progn
 768                            (declaim (inline ,@funs))
 769                            ,@defines))))
 770
 771 (define-twoarg-mpfr-int
 772     (mpfr_sin_cos
 773      mpfr_atan2
 774      mpfr_sinh_cosh
 775      mpfr_agm
 776      mpfr_hypot))
 777
 778 (declaim (inline mpfr_fac_ui
 779                  mpfr_zeta_ui
 780                  mpfr_jn
 781                  mpfr_yn
 782                  mpfr_fma
 783                  mpfr_fms))
 784
 785 (define-alien-routine mpfr_fac_ui int
 786   (r (* (struct mpfrfloat)))
 787   (op unsigned-long)
 788   (rnd mpfr_rnd_enum))
 789
 790 (define-alien-routine mpfr_zeta_ui int
 791   (r (* (struct mpfrfloat)))
 792   (op unsigned-long)
 793   (rnd mpfr_rnd_enum))
 794
 795 (define-alien-routine mpfr_jn int
 796   (r (* (struct mpfrfloat)))
 797   (n long)
 798   (op (* (struct mpfrfloat)))
 799   (rnd mpfr_rnd_enum))
 800
 801 (define-alien-routine mpfr_yn int
 802   (r (* (struct mpfrfloat)))
 803   (n long)
 804   (op (* (struct mpfrfloat)))
 805   (rnd mpfr_rnd_enum))
 806
 807 (define-alien-routine mpfr_fma int
 808   (r (* (struct mpfrfloat)))
 809   (op1 (* (struct mpfrfloat)))
 810   (op2 (* (struct mpfrfloat)))
 811   (op3 (* (struct mpfrfloat)))
 812   (rnd mpfr_rnd_enum))
 813
 814 (define-alien-routine mpfr_fms int
 815   (r (* (struct mpfrfloat)))
 816   (op1 (* (struct mpfrfloat)))
 817   (op2 (* (struct mpfrfloat)))
 818   (op3 (* (struct mpfrfloat)))
 819   (rnd mpfr_rnd_enum))
 820
 821
 822 ;;; constant initialization
 823
 824 (defmacro define-const-mpfr-int (funs)
 825   (loop for i in funs collect `(define-alien-routine ,i int
 826                                  (r (* (struct mpfrfloat)))
 827                                  (rnd mpfr_rnd_enum))
 828           into defines
 829         finally (return `(progn
 830                            (declaim (inline ,@funs))
 831                            ,@defines))))
 832
 833 (define-const-mpfr-int
 834     (mpfr_const_log2
 835      mpfr_const_pi
 836      mpfr_const_euler
 837      mpfr_const_catalan))
 838
 839 (declaim (inline mpfr_sum))
 840 (define-alien-routine mpfr_sum int
 841   (r (* (struct mpfrfloat)))
 842   (tab (* (* (struct mpfrfloat))))
 843   (n unsigned-long)
 844   (rnd mpfr_rnd_enum))
 845
 846
 847 ;;; comparison functions
 848
 849 (declaim (inline mpfr_cmp
 850                  mpfr_cmp_ui
 851                  mpfr_cmp_si
 852                  mpfr_cmp_d
 853                  mpfr_cmp_z
 854                  mpfr_cmp_q
 855                  mpfr_cmp_ui_2exp
 856                  mpfr_cmp_si_2exp
 857                  mpfr_cmpabs))
 858
 859 (define-alien-routine mpfr_cmp int
 860   (op1 (* (struct mpfrfloat)))
 861   (op2 (* (struct mpfrfloat))))
 862
 863 (define-alien-routine mpfr_cmp_ui int
 864   (op1 (* (struct mpfrfloat)))
 865   (op2 unsigned-long))
 866
 867 (define-alien-routine mpfr_cmp_si int
 868   (op1 (* (struct mpfrfloat)))
 869   (op2 long))
 870
 871 (define-alien-routine mpfr_cmp_d int
 872   (op1 (* (struct mpfrfloat)))
 873   (op2 double))
 874
 875 (define-alien-routine mpfr_cmp_z int
 876   (op1 (* (struct mpfrfloat)))
 877   (op2 (* (struct sb-gmp::gmpint))))
 878
 879 (define-alien-routine mpfr_cmp_q int
 880   (op1 (* (struct mpfrfloat)))
 881   (op2 (* (struct sb-gmp::gmprat))))
 882
 883 (define-alien-routine mpfr_cmp_ui_2exp int
 884   (op1 (* (struct mpfrfloat)))
 885   (op2 unsigned-long)
 886   (exp long))
 887
 888 (define-alien-routine mpfr_cmp_si_2exp int
 889   (op1 (* (struct mpfrfloat)))
 890   (op2 long)
 891   (exp long))
 892
 893 (define-alien-routine mpfr_cmpabs int
 894   (op1 (* (struct mpfrfloat)))
 895   (op2 (* (struct mpfrfloat))))
 896
 897
 898 (defmacro define-onearg-mpfr-bool (funs)
 899   (loop for i in funs collect `(define-alien-routine ,i boolean
 900                                  (op (* (struct mpfrfloat))))
 901           into defines
 902         finally (return `(progn
 903                            (declaim (inline ,@funs))
 904                            ,@defines))))
 905
 906 (define-onearg-mpfr-bool
 907     (mpfr_nan_p
 908      mpfr_inf_p
 909      mpfr_number_p
 910      mpfr_zero_p
 911      mpfr_integer_p
 912      mpfr_regular_p))
 913
 914 (declaim (inline mpfr_sgn))
 915 (define-alien-routine mpfr_sgn int
 916   (op (* (struct mpfrfloat))))
 917
 918 (defmacro define-twoarg-mpfr-bool (funs)
 919   (loop for i in funs collect `(define-alien-routine ,i boolean
 920                                  (op1 (* (struct mpfrfloat)))
 921                                  (op2 (* (struct mpfrfloat))))
 922           into defines
 923         finally (return `(progn
 924                            (declaim (inline ,@funs))
 925                            ,@defines))))
 926
 927 (define-twoarg-mpfr-bool
 928     (mpfr_greater_p
 929      mpfr_greaterequal_p
 930      mpfr_less_p
 931      mpfr_lessequal_p
 932      mpfr_equal_p
 933      mpfr_lessgreater_p
 934      mpfr_unordered_p))
 935
 936
 937 ;;; miscellaneous functions
 938
 939 (defmacro define-mpfr-void (funs)
 940   (loop for i in funs collect `(define-alien-routine ,i void)
 941           into defines
 942         finally (return `(progn
 943                            (declaim (inline ,@funs))
 944                            ,@defines))))
 945
 946 (define-mpfr-void
 947     (mpfr_clear_underflow
 948      mpfr_clear_overflow
 949      mpfr_clear_divby0
 950      mpfr_clear_nanflag
 951      mpfr_clear_inexflag
 952      mpfr_clear_erangeflag
 953      mpfr_set_underflow
 954      mpfr_set_overflow
 955      mpfr_set_divby0
 956      mpfr_set_nanflag
 957      mpfr_set_inexflag
 958      mpfr_set_erangeflag
 959      mpfr_clear_flags))
 960
 961 (defmacro define-mpfr-bool (funs)
 962   (loop for i in funs collect `(define-alien-routine ,i boolean)
 963           into defines
 964         finally (return `(progn
 965                            (declaim (inline ,@funs))
 966                            ,@defines))))
 967
 968 (define-mpfr-bool
 969     (mpfr_underflow_p
 970      mpfr_overflow_p
 971      mpfr_divby0_p
 972      mpfr_nanflag_p
 973      mpfr_inexflag_p
 974      mpfr_erangeflag_p))
 975
 976
 977 (declaim (inline mpfr_urandomb
 978                  mpfr_urandom))
 979
 980 (define-alien-routine mpfr_urandomb int
 981   (op (* (struct mpfrfloat)))
 982   (s (* (struct sb-gmp::gmprandstate))))
 983
 984 (define-alien-routine mpfr_urandom int
 985   (op (* (struct mpfrfloat)))
 986   (s (* (struct sb-gmp::gmprandstate)))
 987   (rnd mpfr_rnd_enum))
 988
 989
 990 ;;; integer and remainder related functions / rounding
 991
 992 (define-onearg-mpfr-int
 993     (mpfr_rint
 994      mpfr_rint_ceil
 995      mpfr_rint_floor
 996      mpfr_rint_round
 997      mpfr_rint_trunc
 998      mpfr_frac))
 999
1000 (declaim (inline mpfr_ceil
1001                  mpfr_floor
1002                  mpfr_round
1003                  mpfr_trunc
1004                  mpfr_modf
1005                  mpfr_fmod
1006                  mpfr_remainder
1007                  mpfr_remquo))
1008
1009 (define-alien-routine mpfr_ceil int
1010   (r (* (struct mpfrfloat)))
1011   (op (* (struct mpfrfloat))))
1012
1013 (define-alien-routine mpfr_floor int
1014   (r (* (struct mpfrfloat)))
1015   (op (* (struct mpfrfloat))))
1016
1017 (define-alien-routine mpfr_round int
1018   (r (* (struct mpfrfloat)))
1019   (op (* (struct mpfrfloat))))
1020
1021 (define-alien-routine mpfr_trunc int
1022   (r (* (struct mpfrfloat)))
1023   (op (* (struct mpfrfloat))))
1024
1025 (define-alien-routine mpfr_modf int
1026   (ir (* (struct mpfrfloat)))
1027   (fr (* (struct mpfrfloat)))
1028   (op (* (struct mpfrfloat)))
1029   (rnd mpfr_rnd_enum))
1030
1031 (define-alien-routine mpfr_fmod int
1032   (r (* (struct mpfrfloat)))
1033   (x (* (struct mpfrfloat)))
1034   (y (* (struct mpfrfloat)))
1035   (rnd mpfr_rnd_enum))
1036
1037 (define-alien-routine mpfr_remainder int
1038   (r (* (struct mpfrfloat)))
1039   (x (* (struct mpfrfloat)))
1040   (y (* (struct mpfrfloat)))
1041   (rnd mpfr_rnd_enum))
1042
1043 (define-alien-routine mpfr_remquo int
1044   (r (* (struct mpfrfloat)))
1045   (q (* long))
1046   (x (* (struct mpfrfloat)))
1047   (y (* (struct mpfrfloat)))
1048   (rnd mpfr_rnd_enum))
1049
1050
1051 ;;;; lisp interface
1052
1053 (defparameter *mpfr-rnd* :mpfr_rndn)
1054 (defparameter *mpfr-base* 10)
1055
1056 (declaim (inline mpfr-float-ref
1057                  make-mpfr-float))
1058
1059 (defstruct (mpfr-float (:constructor %make-mpfr-float))
1060   (ref (make-alien (struct mpfrfloat))
1061    :type (alien (* (struct mpfrfloat))) :read-only t))
1062
1063 (defun make-mpfr-float (&optional (precision +mpfr-precision+))
1064   (declare (optimize (speed 3) (space 3)))
1065   (let* ((float (%make-mpfr-float))
1066          (ref (mpfr-float-ref float)))
1067     (mpfr_init2 ref precision)
1068     (sb-ext:finalize float (lambda ()
1069                              (declare (optimize (speed 3) (space 3) (safety 0)))
1070                              (mpfr_clear ref)
1071                              (free-alien ref)))
1072     float))
1073
1074 (defmacro with-precision (value &body body)
1075   `(let* ((old +mpfr-precision+)
1076           (+mpfr-precision+ ,value))
1077      (unwind-protect
1078           (progn
1079             (mpfr_set_default_prec ,value)
1080             ,@body)
1081        (mpfr_set_default_prec old))))
1082
1083 (defun set-precision (value)
1084   (check-type value (unsigned-byte #.sb-vm:n-word-bits))
1085   (mpfr_set_default_prec value)
1086   (setf +mpfr-precision+ value))
1087
1088 ;;; printing and reader syntax
1089
1090 (defmethod print-object ((obj mpfr-float) stream)
1091   (multiple-value-bind (str exp sign)
1092       (mpfr-float-to-string obj)
1093     (declare (type (integer -1 1) sign))
1094     (if *print-readably*
1095         (if (minusp exp)
1096             (format stream "#M~s" str)
1097             (case sign
1098               (0 (format stream "#M\"0\""))
1099               (1 (format stream "#M\"~a@~s\"" str
1100                          (- exp (length str))))
1101               (-1 (format stream "#M\"~a@~s\"" str
1102                           (- (1+ exp) (length str))))))
1103         (if (minusp exp)
1104             (format stream "~a" str)
1105             (case sign
1106               (0 (format stream "0"))
1107               (1 (format stream ".~ae+~s" str exp))
1108               (-1 (format stream "-.~ae+~s" (subseq str 1) exp)))))))
1109
1110 (defun mpfr-float-to-string (x &optional (rnd *mpfr-rnd*))
1111   (let ((xr (mpfr-float-ref x)))
1112     (with-alien ((exp long)
1113                  (str (* char)))
1114       (setf exp -1)
1115       (setf str (mpfr_get_str NIL (addr exp) *print-base* 0 xr rnd))
1116       (multiple-value-prog1
1117           (values (cast str c-string) exp (mpfr_cmp_ui xr 0))
1118         (mpfr_free_str str)))))
1119
1120 (defun mpfr-reader (stream subchar arg)
1121   (declare (ignore subchar arg))
1122   (let ((mpfr (make-mpfr-float)))
1123     (mpfr_set_str (mpfr-float-ref mpfr)
1124                   (read stream t nil t)
1125                   *mpfr-base* *mpfr-rnd*)
1126     mpfr))
1127
1128 (defun enable-mpfr-syntax (readtable)
1129   (set-dispatch-macro-character #\# #\M #'mpfr-reader readtable))
1130 (enable-mpfr-syntax *readtable*)
1131
1132
1133 ;;; arithmetic functions
1134
1135 (defun add (x y &optional (round *mpfr-rnd*))
1136   (if (typep x 'mpfr-float)
1137       (let* ((res (make-mpfr-float))
1138              (r (mpfr-float-ref res))
1139              (xr (mpfr-float-ref x))
1140              (i (etypecase y
1141                   (mpfr-float
1142                    (mpfr_add r xr (mpfr-float-ref y) round))
1143                   ((unsigned-byte #.sb-vm:n-word-bits)
1144                    (mpfr_add_ui r xr y round))
1145                   ((signed-byte #.sb-vm:n-word-bits)
1146                    (mpfr_add_si r xr y round))
1147                   (double-float
1148                    (mpfr_add_d r xr y round))
1149                   (integer
1150                    (sb-gmp::with-mpz-vars ((y gy))
1151                      (mpfr_add_z r xr (addr gy) round)))
1152                   (rational
1153                    (sb-gmp::with-mpq-var (y qy)
1154                      (mpfr_add_q r xr (addr qy) round))))))
1155         (values res i))
1156       (etypecase y
1157         (mpfr-float
1158          (add y x)))))
1159
1160 (defun sub (x y &optional (round *mpfr-rnd*))
1161   (let* ((res (make-mpfr-float))
1162          (r (mpfr-float-ref res))
1163          (i (etypecase x
1164               (mpfr-float
1165                (let ((xr (mpfr-float-ref x)))
1166                  (etypecase y
1167                    (mpfr-float
1168                     (mpfr_sub r xr (mpfr-float-ref y) round))
1169                    ((unsigned-byte #.sb-vm:n-word-bits)
1170                     (mpfr_sub_ui r xr y round))
1171                    ((signed-byte #.sb-vm:n-word-bits)
1172                     (mpfr_sub_si r xr y round))
1173                    (double-float
1174                     (mpfr_sub_d r xr y round))
1175                    (integer
1176                     (sb-gmp::with-mpz-vars ((y gy))
1177                       (mpfr_sub_z r xr (addr gy) round)))
1178                    (rational
1179                     (sb-gmp::with-mpq-var (y qy)
1180                       (mpfr_sub_q r xr (addr qy) round))))))
1181               ((unsigned-byte #.sb-vm:n-word-bits)
1182                (etypecase y
1183                  (mpfr-float
1184                   (mpfr_ui_sub r x (mpfr-float-ref y) round))))
1185               ((signed-byte #.sb-vm:n-word-bits)
1186                (etypecase y
1187                  (mpfr-float
1188                   (mpfr_si_sub r x (mpfr-float-ref y) round))))
1189               (double-float
1190                (etypecase y
1191                  (mpfr-float
1192                   (mpfr_d_sub r x (mpfr-float-ref y) round))))
1193               (integer
1194                (etypecase y
1195                  (mpfr-float
1196                   (sb-gmp::with-mpz-vars ((x gx))
1197                     (mpfr_z_sub r (addr gx) (mpfr-float-ref y) round))))))))
1198     (values res i)))
1199
1200 (defun mul (x y &optional (round *mpfr-rnd*))
1201   (if (typep x 'mpfr-float)
1202       (let* ((res (make-mpfr-float))
1203              (r (mpfr-float-ref res))
1204              (xr (mpfr-float-ref x))
1205              (i (etypecase y
1206                   (mpfr-float
1207                    (mpfr_mul r xr (mpfr-float-ref y) round))
1208                   ((unsigned-byte #.sb-vm:n-word-bits)
1209                    (mpfr_mul_ui r xr y round))
1210                   ((signed-byte #.sb-vm:n-word-bits)
1211                    (mpfr_mul_si r xr y round))
1212                   (double-float
1213                    (mpfr_mul_d r xr y round))
1214                   (integer
1215                    (sb-gmp::with-mpz-vars ((y gy))
1216                      (mpfr_mul_z r xr (addr gy) round)))
1217                   (rational
1218                    (sb-gmp::with-mpq-var (y qy)
1219                      (mpfr_mul_q r xr (addr qy) round))))))
1220         (values res i))
1221       (etypecase y
1222         (mpfr-float
1223          (mul y x)))))
1224
1225 (defun sqare (x &optional (round *mpfr-rnd*))
1226   (let ((r (make-mpfr-float)))
1227     (values r (mpfr_sqr (mpfr-float-ref r) (mpfr-float-ref x) round))))
1228
1229 (defun div (x y &optional (round *mpfr-rnd*))
1230   (let* ((res (make-mpfr-float))
1231          (r (mpfr-float-ref res))
1232          (i (etypecase x
1233               (mpfr-float
1234                (let ((xr (mpfr-float-ref x)))
1235                  (etypecase y
1236                    (mpfr-float
1237                     (mpfr_div r xr (mpfr-float-ref y) round))
1238                    ((unsigned-byte #.sb-vm:n-word-bits)
1239                     (mpfr_div_ui r xr y round))
1240                    ((signed-byte #.sb-vm:n-word-bits)
1241                     (mpfr_div_si r xr y round))
1242                    (double-float
1243                     (mpfr_div_d r xr y round))
1244                    (integer
1245                     (sb-gmp::with-mpz-vars ((y gy))
1246                       (mpfr_div_z r xr (addr gy) round)))
1247                    (rational
1248                     (sb-gmp::with-mpq-var (y qy)
1249                       (mpfr_div_q r xr (addr qy) round))))))
1250               ((unsigned-byte #.sb-vm:n-word-bits)
1251                (etypecase y
1252                  (mpfr-float
1253                   (mpfr_ui_div r x (mpfr-float-ref y) round))))
1254               ((signed-byte #.sb-vm:n-word-bits)
1255                (etypecase y
1256                  (mpfr-float
1257                   (mpfr_si_div r x (mpfr-float-ref y) round))))
1258               (double-float
1259                (etypecase y
1260                  (mpfr-float
1261                   (mpfr_d_div r x (mpfr-float-ref y) round)))))))
1262     (values res i)))
1263
1264 (defun sqrt (x &optional (round *mpfr-rnd*))
1265   (let* ((res (make-mpfr-float))
1266          (r (mpfr-float-ref res))
1267          (i (etypecase x
1268               ((unsigned-byte #.sb-vm:n-word-bits)
1269                (mpfr_sqrt_ui r x round))
1270               (mpfr-float
1271                (mpfr_sqrt r (mpfr-float-ref x) round)))))
1272     (values res i)))
1273
1274 (defun reciprocal-sqrt (x &optional (round *mpfr-rnd*))
1275   (let* ((res (make-mpfr-float))
1276          (r (mpfr-float-ref res))
1277          (i (mpfr_rec_sqrt r (mpfr-float-ref x) round)))
1278     (values res i)))
1279
1280 (defun cubic-root (x &optional (round *mpfr-rnd*))
1281   (let* ((res (make-mpfr-float))
1282          (r (mpfr-float-ref res))
1283          (i (mpfr_cbrt r (mpfr-float-ref x) round)))
1284     (values res i)))
1285
1286 (defun k-root (x k &optional (round *mpfr-rnd*))
1287   (check-type k (unsigned-byte #.sb-vm:n-word-bits))
1288   (let* ((res (make-mpfr-float))
1289          (r (mpfr-float-ref res))
1290          (i (mpfr_root r (mpfr-float-ref x) k round)))
1291     (values res i)))
1292
1293 (defun power (x y &optional (round *mpfr-rnd*))
1294   (let* ((res (make-mpfr-float))
1295          (r (mpfr-float-ref res))
1296          (i (etypecase x
1297               (mpfr-float
1298                (let ((xr (mpfr-float-ref x)))
1299                  (etypecase y
1300                    (mpfr-float
1301                     (mpfr_pow r xr (mpfr-float-ref y) round))
1302                    ((unsigned-byte #.sb-vm:n-word-bits)
1303                     (mpfr_pow_ui r xr y round))
1304                    ((signed-byte #.sb-vm:n-word-bits)
1305                     (mpfr_pow_si r xr y round))
1306                    (integer
1307                     (sb-gmp::with-mpz-vars ((y gy))
1308                       (mpfr_pow_z r xr (addr gy) round))))))
1309               ((unsigned-byte #.sb-vm:n-word-bits)
1310                (etypecase y
1311                  (mpfr-float
1312                   (mpfr_ui_pow r x (mpfr-float-ref y) round)))))))
1313     (values res i)))
1314
1315 (defun negate (x &optional (round *mpfr-rnd*))
1316   (let* ((res (make-mpfr-float))
1317          (r (mpfr-float-ref res))
1318          (i (mpfr_neg r (mpfr-float-ref x) round)))
1319     (values res i)))
1320
1321 (defun abs (x &optional (round *mpfr-rnd*))
1322   (let* ((res (make-mpfr-float))
1323          (r (mpfr-float-ref res))
1324          (i (mpfr_abs r (mpfr-float-ref x) round)))
1325     (values res i)))
1326
1327 (defun dim (x y &optional (round *mpfr-rnd*))
1328   (let* ((res (make-mpfr-float))
1329          (r (mpfr-float-ref res))
1330          (i (mpfr_dim r (mpfr-float-ref x) (mpfr-float-ref y) round)))
1331     (values res i)))
1332
1333 (defun mul-2-raised (x y &optional (round *mpfr-rnd*))
1334   "Compute X*(2^Y)."
1335   (let* ((res (make-mpfr-float))
1336          (r (mpfr-float-ref res))
1337          (i (etypecase y
1338               ((unsigned-byte #.sb-vm:n-word-bits)
1339                (mpfr_mul_2ui r (mpfr-float-ref x) y round))
1340               ((signed-byte #.sb-vm:n-word-bits)
1341                (mpfr_mul_2si r (mpfr-float-ref x) y round)))))
1342     (values res i)))
1343
1344 (defun div-2-raised (x y &optional (round *mpfr-rnd*))
1345   "Compute X/(2^Y)."
1346   (let* ((res (make-mpfr-float))
1347          (r (mpfr-float-ref res))
1348          (i (etypecase y
1349               ((unsigned-byte #.sb-vm:n-word-bits)
1350                (mpfr_div_2ui r (mpfr-float-ref x) y round))
1351               ((signed-byte #.sb-vm:n-word-bits)
1352                (mpfr_div_2si r (mpfr-float-ref x) y round)))))
1353     (values res i)))
1354
1355 ;;; special functions
1356
1357 (defmacro define-onearg-mpfr-funs (funs)
1358   (loop for (clfun mfun) in funs
1359         collect `(defun ,clfun (x &optional (round *mpfr-rnd*))
1360                    (let* ((result (make-mpfr-float))
1361                           (i (,mfun (mpfr-float-ref result)
1362                                     (mpfr-float-ref x)
1363                                     round)))
1364                      (values result i)))
1365           into defines
1366         finally (return `(progn
1367                            ,@defines))))
1368
1369 (define-onearg-mpfr-funs
1370     ((log mpfr_log)
1371      (log2 mpfr_log2)
1372      (log10 mpfr_log10)
1373      (exp mpfr_exp)
1374      (exp2 mpfr_exp2)
1375      (exp10 mpfr_exp10)
1376      (cos mpfr_cos)
1377      (sin mpfr_sin)
1378      (tan mpfr_tan)
1379      (sec mpfr_sec)
1380      (csc mpfr_csc)
1381      (cot mpfr_cot)
1382      (acos mpfr_acos)
1383      (asin mpfr_asin)
1384      (cosh mpfr_cosh)
1385      (sinh mpfr_sinh)
1386      (tanh mpfr_tanh)
1387      (sech mpfr_sech)
1388      (csch mpfr_csch)
1389      (coth mpfr_coth)
1390      (acosh mpfr_acosh)
1391      (asinh mpfr_asinh)
1392      (atanh mpfr_atanh)
1393      (log1p mpfr_log1p)
1394      (expm1 mpfr_expm1)
1395      (eint mpfr_eint)
1396      (li2 mpfr_li2)
1397      (gamma mpfr_gamma)
1398      (lngamma mpfr_lngamma)
1399      (digamma mpfr_digamma)
1400      (erf mpfr_erf)
1401      (erfc mpfr_erfc)
1402      (j0 mpfr_j0)
1403      (j1 mpfr_j1)
1404      (y0 mpfr_y0)
1405      (y1 mpfr_y1)
1406      (ai mpfr_ai)))
1407
1408 (defun atan (y &optional x (round *mpfr-rnd*))
1409   (if x
1410       (let* ((result (make-mpfr-float))
1411              (i (mpfr_atan2 (mpfr-float-ref result)
1412                             (mpfr-float-ref y)
1413                             (mpfr-float-ref x)
1414                             round)))
1415         (values result i))
1416       (let* ((result (make-mpfr-float))
1417              (i (mpfr_atan (mpfr-float-ref result)
1418                            (mpfr-float-ref y)
1419                            round)))
1420         (values result i))))
1421
1422 (defun sin-cos (x &optional (round *mpfr-rnd*))
1423   (let* ((sin (make-mpfr-float))
1424          (cos (make-mpfr-float))
1425          (i (mpfr_sin_cos (mpfr-float-ref sin)
1426                           (mpfr-float-ref cos)
1427                           (mpfr-float-ref x)
1428                           round)))
1429     (values sin cos i)))
1430
1431 (defun sinh-cosh (x &optional (round *mpfr-rnd*))
1432   (let* ((sin (make-mpfr-float))
1433          (cos (make-mpfr-float))
1434          (i (mpfr_sinh_cosh (mpfr-float-ref sin)
1435                             (mpfr-float-ref cos)
1436                             (mpfr-float-ref x)
1437                             round)))
1438     (values sin cos i)))
1439
1440 (defun arithmetic-geometric-mean (u v &optional (round *mpfr-rnd*))
1441   (let* ((result (make-mpfr-float))
1442          (i (mpfr_agm (mpfr-float-ref result)
1443                       (mpfr-float-ref u)
1444                       (mpfr-float-ref v)
1445                       round)))
1446     (values result i)))
1447
1448 (defun hypot (x y &optional (round *mpfr-rnd*))
1449   (let* ((result (make-mpfr-float))
1450          (i (mpfr_hypot (mpfr-float-ref result)
1451                         (mpfr-float-ref x)
1452                         (mpfr-float-ref y)
1453                         round)))
1454     (values result i)))
1455
1456 (defun fac (x &optional (round *mpfr-rnd*))
1457   (let* ((result (make-mpfr-float))
1458          (i (mpfr_fac_ui (mpfr-float-ref result)
1459                          x
1460                          round)))
1461     (values result i)))
1462
1463 (defun zeta (x &optional (round *mpfr-rnd*))
1464   (let* ((result (make-mpfr-float))
1465          (i (etypecase x
1466               (mpfr-float
1467                (mpfr_zeta (mpfr-float-ref result)
1468                           (mpfr-float-ref x)
1469                           round))
1470               ((unsigned-byte #.sb-vm:n-word-bits)
1471                (mpfr_zeta_ui (mpfr-float-ref result)
1472                              x
1473                              round)))))
1474     (values result i)))
1475
1476 (defun jn (n x &optional (round *mpfr-rnd*))
1477   (check-type n (unsigned-byte #.sb-vm:n-word-bits))
1478   (let* ((result (make-mpfr-float))
1479          (i (mpfr_jn (mpfr-float-ref result)
1480                      n
1481                      (mpfr-float-ref x)
1482                      round)))
1483     (values result i)))
1484
1485 (defun yn (n x &optional (round *mpfr-rnd*))
1486   (check-type n (unsigned-byte #.sb-vm:n-word-bits))
1487   (let* ((result (make-mpfr-float))
1488          (i (mpfr_yn (mpfr-float-ref result)
1489                      n
1490                      (mpfr-float-ref x)
1491                      round)))
1492     (values result i)))
1493
1494 (defun fma (x y z &optional (round *mpfr-rnd*))
1495   "fma X Y Z = (X * Y) + Z"
1496   (let* ((result (make-mpfr-float))
1497          (i (mpfr_fma (mpfr-float-ref result)
1498                       (mpfr-float-ref x)
1499                       (mpfr-float-ref y)
1500                       (mpfr-float-ref z)
1501                       round)))
1502     (values result i)))
1503
1504 (defun fms (x y z &optional (round *mpfr-rnd*))
1505   "fma X Y Z = (X * Y) - Z"
1506   (let* ((result (make-mpfr-float))
1507          (i (mpfr_fms (mpfr-float-ref result)
1508                       (mpfr-float-ref x)
1509                       (mpfr-float-ref y)
1510                       (mpfr-float-ref z)
1511                       round)))
1512     (values result i)))
1513
1514 (defun sum (seq &optional (round *mpfr-rnd*))
1515   (let ((length (length seq))
1516         (idx -1))
1517     (declare (type (integer -1 #.most-positive-fixnum) idx))
1518     (let ((result (make-mpfr-float))
1519           (ar (make-alien (* (struct mpfrfloat)) length)))
1520       (unwind-protect
1521            (progn
1522              (map nil (lambda (x)
1523                         (setf (deref ar (incf idx))
1524                               (mpfr-float-ref x)))
1525                   seq)
1526              (let ((i (mpfr_sum (mpfr-float-ref result)
1527                                 ar
1528                                 length
1529                                 round)))
1530                (values result i)))
1531         (free-alien ar)))))
1532
1533
1534 ;;; constant values
1535
1536 (defmacro define-const-mpfr-funs (funs)
1537   (loop for (fname mname) in funs
1538         collect `(defun ,fname (&optional (round *mpfr-rnd*))
1539                    (let* ((result (make-mpfr-float))
1540                           (i (,mname (mpfr-float-ref result)
1541                                      round)))
1542                      (values result i)))
1543           into defines
1544         finally (return `(progn
1545                            ,@defines))))
1546
1547 (define-const-mpfr-funs
1548     ((const-log2 mpfr_const_log2)
1549      (const-pi mpfr_const_pi)
1550      (const-euler mpfr_const_euler)
1551      (const-catalan mpfr_const_catalan)))
1552
1553 ;;; comparison functions and predicates
1554
1555 (defmacro define-onearg-mpfr-predicates (funs)
1556   (loop for (fname mname) in funs
1557         collect `(defun ,fname (x)
1558                    (,mname (mpfr-float-ref x)))
1559           into defines
1560         finally (return `(progn
1561                            ,@defines))))
1562
1563 (define-onearg-mpfr-predicates
1564     ((nan-p mpfr_nan_p)
1565      (infinityp mpfr_inf_p)
1566      (numberp mpfr_number_p)
1567      (zerop mpfr_zero_p)
1568      (integerp mpfr_integer_p)
1569      (regularp mpfr_regular_p)))
1570
1571 (defmacro define-twoarg-mpfr-predicates (funs)
1572   (loop for (fname mname) in funs
1573         collect `(defun ,fname (x y)
1574                    (,mname (mpfr-float-ref x)
1575                            (mpfr-float-ref y)))
1576           into defines
1577         finally (return `(progn
1578                            ,@defines))))
1579
1580 (define-twoarg-mpfr-predicates
1581     ((compare-abs mpfr_cmpabs)
1582      (> mpfr_greater_p)
1583      (>= mpfr_greaterequal_p)
1584      (< mpfr_less_p)
1585      (<= mpfr_lessequal_p)
1586      (= mpfr_equal_p)
1587      (/= mpfr_lessgreater_p)
1588      (unorderedp mpfr_unordered_p)))
1589
1590 (defun compare (x y)
1591   (if (typep x 'mpfr-float)
1592       (etypecase y
1593         (mpfr-float
1594          (mpfr_cmp (mpfr-float-ref x)
1595                    (mpfr-float-ref y)))
1596         ((unsigned-byte #.sb-vm:n-word-bits)
1597          (mpfr_cmp_ui (mpfr-float-ref x) y))
1598         ((signed-byte #.sb-vm:n-word-bits)
1599          (mpfr_cmp_si (mpfr-float-ref x) y))
1600         (double-float
1601          (mpfr_cmp_d (mpfr-float-ref x) y))
1602         (integer
1603          (sb-gmp::with-mpz-vars ((y gy))
1604            (mpfr_cmp_z (mpfr-float-ref x) (addr gy))))
1605         (rational
1606          (sb-gmp::with-mpq-var (y qy)
1607            (mpfr_cmp_q (mpfr-float-ref x) (addr qy)))))
1608       (etypecase y
1609         (mpfr-float
1610          (compare y x)))))
1611
1612 (defun compare-2exp (x y exp)
1613   (if (typep x 'mpfr-float)
1614       (etypecase y
1615         ((unsigned-byte #.sb-vm:n-word-bits)
1616          (mpfr_cmp_ui_2exp (mpfr-float-ref x) y exp))
1617         ((signed-byte #.sb-vm:n-word-bits)
1618          (mpfr_cmp_si_2exp (mpfr-float-ref x) y exp)))
1619       (etypecase y
1620         (mpfr-float
1621          (compare-2exp y x exp)))))
1622
1623
1624 ;;; miscellaneous functions
1625
1626 (defmacro define-mpfr-flag-funs (funs)
1627   (loop for (pname mname) in funs
1628         collect `(defun ,pname ()
1629                    (declare (optimize (speed 3) (space 3)))
1630                    (,mname) (values))
1631           into defines
1632         collect pname into names
1633         finally (return `(progn
1634                            (declaim (inline ,@names))
1635                            ,@defines))))
1636
1637 (define-mpfr-flag-funs
1638     ((clear-underflow mpfr_clear_underflow)
1639      (clear-overflow mpfr_clear_overflow)
1640      (clear-div-by-zero mpfr_clear_divby0)
1641      (clear-nan-flag mpfr_clear_nanflag)
1642      (clear-inex-flag mpfr_clear_inexflag)
1643      (clear-erange-flag mpfr_clear_erangeflag)
1644      (set-underflow-flag mpfr_set_underflow)
1645      (set-overflow-flag mpfr_set_overflow)
1646      (set-div-by-zero-flag mpfr_set_divby0)
1647      (set-nan-flag mpfr_set_nanflag)
1648      (set-inex-flag mpfr_set_inexflag)
1649      (set-erange-flag mpfr_set_erangeflag)
1650      (clear-flags mpfr_clear_flags)))
1651
1652
1653 (defmacro define-mpfr-flag-predicates (funs)
1654   (loop for (pname mname) in funs
1655         collect `(defun ,pname ()
1656                    (declare (optimize (speed 3) (space 3)))
1657                    (,mname))
1658           into defines
1659         collect pname into names
1660         finally (return `(progn
1661                            (declaim (inline ,@names))
1662                            ,@defines))))
1663
1664 (define-mpfr-flag-predicates
1665     ((underflowp mpfr_underflow_p)
1666      (overflowp mpfr_overflow_p)
1667      (div-by-zero-p mpfr_divby0_p)
1668      (nanflag-p mpfr_nanflag_p)
1669      (inexflag-p mpfr_inexflag_p)
1670      (erangeflag-p mpfr_erangeflag_p)))
1671
1672
1673 ;;; random number generation
1674
1675 (defun urandomb (state)
1676   (check-type state sb-gmp::gmp-rstate)
1677   (let* ((ref (sb-gmp::gmp-rstate-ref state))
1678          (result (make-mpfr-float))
1679          (i (mpfr_urandomb (mpfr-float-ref result)
1680                            ref)))
1681     (values result i)))
1682
1683 (defun urandom (state &optional (round *mpfr-rnd*))
1684   (check-type state sb-gmp::gmp-rstate)
1685   (let* ((ref (sb-gmp::gmp-rstate-ref state))
1686          (result (make-mpfr-float))
1687          (i (mpfr_urandom (mpfr-float-ref result)
1688                           ref
1689                           round)))
1690     (values result i)))
1691
1692
1693 ;;; integer and remainder related functions / rounding
1694
1695 (define-onearg-mpfr-funs
1696     ((rounded-int mpfr_rint)
1697      (rounded-int-ceiling mpfr_rint_ceil)
1698      (rounded-int-floor mpfr_rint_floor)
1699      (rounded-int-round mpfr_rint_round)
1700      (rounded-int-truncate mpfr_rint_trunc)
1701      (fractional mpfr_frac)))
1702
1703 (defmacro define-onearg-no-rnd-mpfr-funs (funs)
1704   (loop for (clfun mfun) in funs
1705         collect `(defun ,clfun (x)
1706                    (let* ((result (make-mpfr-float))
1707                           (i (,mfun (mpfr-float-ref result)
1708                                     (mpfr-float-ref x))))
1709                      (values result i)))
1710           into defines
1711         finally (return `(progn
1712                            ,@defines))))
1713
1714 (define-onearg-no-rnd-mpfr-funs
1715     ((ceil mpfr_ceil)
1716      (floor mpfr_floor)
1717      (round mpfr_round)
1718      (truncate mpfr_trunc)))
1719
1720 (defun modf (x &optional (round *mpfr-rnd*))
1721   (let* ((integral (make-mpfr-float))
1722          (fractional (make-mpfr-float))
1723          (i (mpfr_modf (mpfr-float-ref integral)
1724                        (mpfr-float-ref fractional)
1725                        (mpfr-float-ref x)
1726                        round)))
1727     (values integral fractional i)))
1728
1729 (defun fmod (x y &optional (round *mpfr-rnd*))
1730   (let* ((result (make-mpfr-float))
1731          (i (mpfr_fmod (mpfr-float-ref result)
1732                        (mpfr-float-ref x)
1733                        (mpfr-float-ref y)
1734                        round)))
1735     (values result i)))
1736
1737 (defun remainder (x y &optional (round *mpfr-rnd*))
1738   (let* ((result (make-mpfr-float))
1739          (i (mpfr_remainder (mpfr-float-ref result)
1740                             (mpfr-float-ref x)
1741                             (mpfr-float-ref y)
1742                             round)))
1743     (values result i)))
1744
1745 (defun remainder-quot (x y &optional (round *mpfr-rnd*))
1746   (with-alien ((q long))
1747     (let* ((result (make-mpfr-float))
1748            (i (mpfr_remquo (mpfr-float-ref result)
1749                            (addr q)
1750                            (mpfr-float-ref x)
1751                            (mpfr-float-ref y)
1752                            round)))
1753       (values result q i))))
1754
1755
1756 ;;; conversion
1757
1758 (declaim (inline mpfr_fits_ulong_p mpfr_fits_slong_p))
1759
1760 (define-alien-routine mpfr_fits_ulong_p boolean
1761   (x (* (struct mpfrfloat)))
1762   (rnd mpfr_rnd_enum))
1763
1764 (define-alien-routine mpfr_fits_slong_p boolean
1765   (x (* (struct mpfrfloat)))
1766   (rnd mpfr_rnd_enum))
1767
1768 (defun coerce (x type &optional (round *mpfr-rnd*))
1769   (cond
1770     ((typep x 'mpfr-float)
1771      (let ((x-ref (mpfr-float-ref x)))
1772        (case type
1773          (single-float
1774           (mpfr_get_flt x-ref round))
1775          (double-float
1776           (mpfr_get_d x-ref round))
1777          (mpfr-float
1778           (let ((result (make-mpfr-float)))
1779             (mpfr_set (mpfr-float-ref result) x-ref round)
1780             result))
1781          (integer
1782           (unless (numberp x)
1783             (error "Cannot coerce ~s to ~s. Argument must be an actual number."
1784                    x type))
1785           (cond
1786             ((mpfr_fits_slong_p x-ref round)
1787              (mpfr_get_si x-ref round))
1788             ((mpfr_fits_ulong_p x-ref round)
1789              (mpfr_get_ui x-ref round))
1790             (t
1791              (sb-gmp::with-gmp-mpz-results (rop)
1792                (mpfr_get_z (addr rop) x-ref round)))))
1793          (t
1794           (error "TYPE must be one of SINGLE-FLOAT, DOUBLE-FLOAT or INTEGER.")))))
1795     ((eql type 'mpfr-float)
1796      (let ((result (make-mpfr-float)))
1797        (etypecase x
1798          ((signed-byte #.sb-vm:n-word-bits)
1799           (mpfr_set_si (mpfr-float-ref result) x round))
1800          ((unsigned-byte #.sb-vm:n-word-bits)
1801           (mpfr_set_ui (mpfr-float-ref result) x round))
1802          (integer
1803           (sb-gmp::with-mpz-vars ((x gx))
1804             (mpfr_set_z (mpfr-float-ref result) (addr gx) round)))
1805          (single-float
1806           (mpfr_set_flt (mpfr-float-ref result) x round))
1807          (double-float
1808           (mpfr_set_d (mpfr-float-ref result) x round))
1809          (ratio
1810           (sb-gmp::with-mpq-var (x qx)
1811             (mpfr_set_q (mpfr-float-ref result) (addr qx) round))))
1812        result))
1813     (t
1814      (error "Unable to handle the combination of ~S and ~S." x type))))