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