contrib/sb-mpfr/mpfr.lisp

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