(     Title:  FORTH-GMPFR LIBRARY
       File:  gmpfr.fs
     Author:  David N. Williams
    Version:  0.9.2
    License:  LGPL 3
  Test file:  gmpfr-test.fs
   Log file:  gmpfr.log
    Revised:  February 6, 2021

All code derived from elsewhere is believed to be either compatible with
the LGPL, or in the public domain.  Code neither derived from elsewhere
nor in the public domain is
)
\ Copyright (C) 2011, 2012, 2018-2021 David N. Williams
(
This library is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 3 of the License, or at
your option any later version.

This library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with this library; if not, write to the Free Software Foundation,
Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA.

If you take advantage of the option in the LGPL to put a particular
version of this library under the GPL, the author[s] would regard it as
polite if you would put any direct modifications under the LGPL as well,
and include a copy of this request near the beginning of the modified
library source.  A "direct modification" is one that enhances or extends
the library in line with its original concept, as opposed to developing
a distinct application or library which might use it.

This file defines Forth words for gmp big integers, rationals, and
floats, and mpfr floats, designated b, q, mf, and rf for short, in terms
of bindings for the gmp library, libgmp 6.2.1, the mpfr library, libmpfr
4.1.0, and the garbage-collected record stack library, gcrstack.fs.

Currently the documentation of word actions is skimpy or nonexistent. 
It had been hoped that the rather well written GMP library documentation
could serve that purpose, but now we know better.  It is nevertheless
indispensable:

  https://gmplib.org/gmp-man-6.2.1.pdf
  https://www.mpfr.org/mpfr-current/mpfr.pdf

The test file gmpfr-test.fs contains many syntax examples.  It tests the
binding scheme, and is not a comprehensive test of the correctness of
the GNU libraries themselves.

Garbage collection is supposed to be transparent to the user. See,
however, the warning in gcrstack.fs about words that fetch references
from the record stack to the data stack.  It applies to the stack fetch
words in this library.  Generally speaking, such references should
appear on the data stack only under the hood.

The current design puts a further restriction on the use of such words,
due to the immediate freeing of the memory managed by gmp and mpfr
records when they are garbaged.  Operations that could result in
garbaging, such as drops, must be postponed until after gmp functions
that use the record data are finished. Again, such operations are
expected to occur only under the hood.

This file provides five independent word sets, all of which are loaded
by default.  An external program can prevent any subset of the five from
loading by predefining the relevant CAN-PROVIDE items listed below as
false constants.  For example,

  false CONSTANT MULTI-RATIONALS
  false CONSTANT MULTI-FLOATS
  INCLUDE gmpfr.fs

would load only the word sets for MULTI-INTEGERS, RELIABLE-FLOATS, AND
RANDOMS.
)

: CAN-PROVIDE  ( "<flag-constant>" -- )
(
If the constant is already defined, parse it away.  Else define it as
true.  A defining version of Bruce McFarling's [CHECK].
)
  >in @ postpone [UNDEFINED]
  IF >in ! true constant ELSE drop THEN ;

CAN-PROVIDE MULTI-INTEGERS
CAN-PROVIDE MULTI-RATIONALS
CAN-PROVIDE MULTI-FLOATS
CAN-PROVIDE RELIABLE-FLOATS
CAN-PROVIDE RANDOMS

MULTI-INTEGERS MULTI-RATIONALS or MULTI-FLOATS or
RELIABLE-FLOATS or RANDOMS or 0= [IF]
cr .( ***No forth-gmpfr word sets loaded) cr ABORT [THEN] 

s" FORTH-NAME" ENVIRONMENT? dup [IF]
  drop s" pfe" compare 0=
  false REDEFINED-MSG !
[THEN] CONSTANT PFE-HOST

s" gforth" environment? dup [IF]
  ( s flag) nip nip
  0 WARNINGS !
[THEN] CONSTANT GFORTH-HOST

s" IFORTH" environment? dup [IF]
  ( flag) drop
  0 WARNING !
[THEN] CONSTANT IFORTH-HOST

\ libgmp is always loaded
PFE-HOST    [IF] LOADM libgmp-pfe  LOADM zchar  [THEN]
GFORTH-HOST [IF] s" libgmp-gforth.fs"  REQUIRED [THEN]
IFORTH-HOST [IF] s" libgmp-iforth.frt" INCLUDED [THEN]

IFORTH-HOST [IF]
s" gcrstack.fs" INCLUDED
s" mstrings.fs" INCLUDED  \ loads parsing.fs
[ELSE]
s" gcrstack.fs" REQUIRED
s" mstrings.fs" REQUIRED  \ loads parsing.fs
[THEN]

: FORTH-GMPFR-VERSION-S  s" 0.9.2" ;


\ *** TABLE OF CONTENTS
---
  0 NOTATION AND TERMINOLOGY

  1 GENERAL SUPPORT
      MZBUF  /MZBUF  MZ"  -mz-ws
      DEFAULT-MZBUF  DEFAULT-/MZBUF

  2 INTEGERS

  2.1 RECORD SPACE
      DEFAULT-/BSPACE  BSPACE
      make-bspace  clear-bspace  free-bspace

  2.2 GARBAGE COLLECTION
      collect-bgarbage  bgarbage?
      bgc-off  bgc-on  bgc-lock@  bgc-lock!

  2.3 STACK
      (b:  bunused  bdepth  bs@  bs2@  bs3@  >bs  >bs-copy
      BVARIABLE  BCONSTANT  free-bconstants
      b"  b.  b@  b!  b,  bdrop  b2drop  b3drop  bdup  b2dup  bswap
      bnip  bover  btuck  brot  -brot  bpick  bexchange

  2.4 ACCUMULATOR ARITHMETIC
      b>ba  ba>b  by>b   bz>b
      0>ba  0>by  0>bz   n>ba  u>ba  df>ba
      ba>n  ba>u  ba>df  ba>df-exp  ba.

      \ elementary
      -ba  |ba|
      ba+  ba-   bau+  bau-  uba-
      ba*  ban*  bau*  b*a+  bu*a+  ba*-  bu*a-
      ba2^u*

      \ ceiled division
      bac/   bacrem   bac/rem
      bacu/  bacurem  bacu/rem  bac/urem  ba2^uc/  ba2^ucrem

      \ floored division
      baf/   bafrem   baf/rem
      bafu/  bafurem  bafu/rem  baf/urem  ba2^uf/  ba2^ufrem

      \ symmetric division
      bas/   basrem   bas/rem
      basu/  basurem  basu/rem  bas/urem  ba2^us/  ba2^usrem

      \ exact division
      ba/exact  bau/exact

      \ modulo
      bamod  baumod

      \ powers
      ba-pow-mod  ba-pow-mod-sec  ba^u-mod  ba^u  u^u>ba

      \ roots
      ba-nroot  ba-nroot-rem  ba-sqrt  ba-sqrt-rem

      \ number theoretics
      ba-next-prime  ba-gcd  bau-gcd  ba-gcd-st  ba-gcd-s
      ba-lcm  bau-lcm  ba-mod-invert  ba-remove
      fac>ba  ba-bin   bin>ba
      fib>ba  fib2>ba  luc>ba  luc2>ba 

      \ logic and bits
      ba-and  ba-or  ba-xor  ba-invert
      ba-bit-on  ba-bit-off  ba-bit-flip  ba-bit-test

  2.5 ARITHMETIC
      0>b  n>b  u>b  df>b
      b>n  b>u  b>df  b>df-exp
      
      \ elementary
      -b  |b|
      b+   b-   bu+   bu-  ub-  b*  bn*  bu*
      b*+  b*-  bu*+  bu*-
      b2^u*

      \ ceiled division
      bc/   bcrem   bc/rem
      bcu/  bcurem  bcu/rem  bc/urem  b2^uc/  b2^ucrem

      \ floored division
      bf/   bfrem   bf/rem
      bfu/  bfurem  bfu/rem  bf/urem  b2^uf/  b2^ufrem

      \ symmetric division
      bs/   bsrem   bs/rem
      bsu/  bsurem  bsu/rem  bs/urem  b2^us/  b2^usrem

      \ exact division
      b/exact  bu/exact

      \ modulo
      bmod  bumod

      \ divisibility and congruence
      b-divisible  bu-divisible  b2^u-divisible
      b-congruent  bu-congruent  b2^u-congruent

      \ powers
      b-pow-mod  b-pow-mod-sec  b^u-mod  b^u  u^u>b

      \ roots
      b-nroot  b-nroot-rem  b-sqrt  b-sqrt-rem
      b-perfect-pow  b-perfect-squ

      \ number theoretics
      b-prob-prime  b-next-prime
      b-gcd  bu-gcd>u  bu-gcd  b-gcd-st  b-gcd-s
      b-lcm  bu-lcm  b-mod-invert
      b-jacobi  b-legendre  b-kronecker
      bn-kronecker  bu-kronecker  nb-kronecker  ub-kronecker
      b-remove  fac>b  bu-bin  bin>b  fib>b  fib2>b  luc>b  luc2>b

      \ logic and bits
      b-and  b-or  b-xor  b-invert
      b-bit-on  bbit-off  b-bit-flip  b-bit-test
      b-pop-count  b-ham-dist  b-scan0  b-scan1

  2.6 COMPARISONS
      bsgn
      b>      b=      b<      b>=      b<= 
      b0>     b0=     b0<     b0>=     b0<=
      bn>     bn=     bn<     bn>=     bn<=
      bu>     bu=     bu<     bu>=     bu<=
      bdf>    bdf=    bdf<    bdf>=    bdf<=
      b|>|    b|=|    b|<|    b|>=|    b|<=|
      bu|>|   bu|=|   bu|<|   bu|>=|   bu|<=| 
      bdf|>|  bdf|=|  bdf|<|  bdf|>=|  bdf|<=|

  3 RATIONALS

  3.1 RECORD SPACE
      DEFAULT-/QSPACE  QSPACE
      make-qspace  clear-qspace  free-qspace

  3.2 GARBAGE COLLECTION
      collect-qgarbage  qgarbage?
      qgc-off  qgc-on  qgc-lock@  qgc-lock!

  3.3 STACK
      (q:  qunused  qdepth  qs@  qs2@  qs3@  >qs  >qs-copy
      QVARIABLE  QCONSTANT  free-qconstants
      q"  q.  q@  q!  q,  qdrop  q2drop  q3drop  qdup  q2dup  qswap
      qnip  qover  qtuck  qrot  -qrot  qpick  qexchange

  3.4 ACCUMULATOR ARITHMETIC
      q>qa  qa>q  0>qa  nu>qa  uu>qa  df>qa  qa>df  qacanon  qa.

      \ elementary
      -qa  |qa|	 1/qa  qa+  qa-	 qa*  qa/  qa2^u*  qa2^u/

  3.5 ARITHMETIC
      0>q  nu>q  uu>q  df>q  q>df  qcanon

      \ elementary
      -q  |q|  1/q  q+  q-  q*  q/  q2^u*  q2^u/

  3.6 COMPARISONS
      qsgn
      q>      q=      q<      q>=      q<= 
      q0>     q0=     q0<     q0>=     q0<=
      qnu>    qnu=    qnu<    qnu>=    qnu<=
      quu>    quu=    quu<    quu>=    quu<=

  4 FLOATING POINT

  4.1 PRECISION
      set-default-mfprecision  get-default-mfprecision
      set-output-mfprecision   get-output-mfprecision

  4.2 RECORD SPACE
      DEFAULT-/MFSPACE  MFSPACE
      make-mfspace  clear-mfspace  free-mfspace
      mfa-start-precision  mfa-init-precision

  4.3 GARBAGE COLLECTION
      collect-mfgarbage  mfgarbage?
      mfgc-off  mfgc-on  mfgc-lock@  mfgc-lock!

  4.4 STACK
      (mf:  mfunused  mfdepth  mfs@  mfs2@  mfs3@  >mfs  >mfs-copy
      MFVARIABLE  MFCONSTANT  free-mfconstants
      mf"  (mf.)  mf.  mf@  mf!  mf,  mfdrop  mf2drop  mf3drop
      mfdup  mf2dup  mfswap  mfnip  mfover  mftuck  mfrot  -mfrot
      mfpick  mfexchange

  4.5 ACCUMULATOR ARITHMETIC
      mf>mfa  mfa>mf  0>mfa
      n>mfa   u>mfa   df>mfa
      mfa>n   mfa>u   mfa>df  mfa>df-exp  mfa.
      set-mfa-precision  get-mfa-precision

      \ elementary
      -mfa  |mfa|  mfa+  mfa-  mfau+  mfau-  umfa-  mfa*  mfau*
      mfa/  mfau/  umfa/
      mfa-sqrt  usqrt>mfa  mfa^u  mfa2^u*  mfa2^u/  mfa-reldiff

      \ integer rounding
      mfa-ceil  mfa-floor  mfa-trunc  mfa-integer

  4.6 ARITHMETIC
      0>mf  n>mf  u>mf  df>mf
      mf>n  mf>u  mf>df  mf>df-exp

      \ elementary
      -mf  |mf|  mf+  mf-  mfu+  mfu-  umf-  mf*  mfu*  mf/  mfu/  umf/
      mf-sqrt  usqrt>mf  mf^u  mf2^u*  mf2^u/  mf-reldiff

      \ integer rounding
      mf-ceil  mf-floor  mf-trunc  mf-integer

  4.7 COMPARISONS
      mfsgn
      mf>    mf=    mf<    mf>=    mf<= 
      mf0>   mf0=   mf0<   mf0>=   mf0<=
      mfn>   mfn=   mfn<   mfn>=   mfn<=
      mfu>   mfu=   mfu<   mfu>=   mfu<=
      mfdf>  mfdf=  mfdf<  mfdf>=  mfdf<=

  5 RELIABLE FLOATING POINT

  5.1 PRECISION AND RANGE
      set-default-rfprecision get-default-rfprecision
      set-output-rfprecision  get-output-rfprecision
      get-rfemin  get-rfemax  set-rfemin  set-rfemax
      RFEMIN-MIN  RFEMIN-MAX  RFEMAX-MIN  RFEMAX-MAX
      set-IEEE-754

  5.2 ROUNDING
      get-rfround  set-rfround
      RF-NEAREST  RF-UP  RF-DOWN  RF-TOWARD  RF-AWAY

  5.3 EXCEPTION FLAGS
      rfunderflow     rfoverflow      rfinvalid	     rfinexact	    rferange
      set-rfunderflow set-rfoverflow  set-rfinvalid  set-rfinexact  set-rferange
      clr-rfunderflow clr-rfoverflow  clr-rfinvalid  clr-rfinexact  clr-rferange
      clr-rfflags

  5.4 RECORD SPACE
      DEFAULT-/RFSPACE  RFSPACE
      make-rfspace  clear-rfspace  free-rfspace

  5.5 GARBAGE COLLECTION
      collect-rfgarbage  rfgarbage?
      rfgc-off  rfgc-on  rfgc-lock@  rfgc-lock!

  5.6 STACK
      (rf:  rfunused  rfdepth  rfs@  rfs2@  rfs3@  >rfs  >rfs-copy
      RFVARIABLE  RFCONSTANT  free-rfconstants
      rf"  (rf.)  rf.  rf@  rf!  rf,  rfdrop  rf2drop  rf3drop
      rfdup  rf2dup  rfswap  rfnip  rfover  rftuck  rfrot  -rfrot
      rfpick  rfexchange

  5.7 ACCUMULATOR ARITHMETIC
      rf>rfa  rfa>rf
      0>rfa   -0>rfa  nan>rfa  -nan>rfa  inf>rfa  -inf>rfa
      n>rfa   u>rfa   df>rfa
      pi>rfa  ln2>rfa  euler>rfa  catalan>rfa 
      rfa>n  rfa>u  rfa>df  rfa>df-exp
      set-rfa-precision  get-rfa-precision  (rfa.)  rfa.

      \ elementary
      -rfa  |rfa|
      rfa+  rfan+  rfau+  rfadf+
      rfa-  rfan-  nrfa-  rfau-   urfa-  rfadf-  dfrfa-
      rfa*  rfan*  rfau*  rfadf*  rfa*+  rfa*-
      rfa/  rfan/  rfau/  nrfa/   urfa/  rfadf/  dfrfa/
      rfa-sqrt  usqrt>rfa  1/rfa-sqrt  rfa-cbrt  rfa-uroot
      rfa-square  rfa-pow  rfa^n  rfa^u  u^rfa  u^u>rfa
      rfa2^n*  rfa2^u*  rfa2^n/  rfa2^u/
      rfa-max  rfa-min  rfa-reldiff  rfa-dim

      \ integer rounding
      rfa-ceil  rfa-floor  rfa-trunc  rfa-integer

      \ special functions
      rfa-ln    rfa-log2  rfa-log10  rfa-exp    rfa-exp2  rfa-exp10
      rfa-cos   rfa-sin   rfa-tan    rfa-sincos rfa-sec   rfa-csc  rfa-cot
      rfa-acos  rfa-asin  rfa-atan   rfa-atan2
      rfa-cosh  rfa-sinh  rfa-tanh   rfa-sinhcosh
      rfa-sech  rfa-csch  rfa-coth   rfa-acosh  rfa-asinh rfa-atanh
      ufac>rfa  rfa-ln1p  rfa-expm1  rfa-eint   rfa-li2
      rfa-gamma rfa-lngamma  rfa-lgamma  rfa-digamma
      rfa-zeta  uzeta>rfa  rfa-erf  rfa-erfc
      rfa-j0  rfa-j1  rfa-jn  rfa-y0  rfa-y1  rfa-yn  rfa-ai
      rfa-agm  rfa-hypot

      \ classification
      rfa-finite  rfa-infinite  rfa-nan  rfa-regular   

      \ data manipulation
      rfa-signbit    rfa-setsign  rfa-copysign
      rfa-nextafter  rfa-nextup   rfa-nextdown
      rfa-checkrange rfa-subnormalize
      
  5.8 ARITHMETIC
      0>rf  -0<rf  nan>rf  -nan>rf  inf>rf  -inf>rf
      n>rf   u>rf  df>rf
      rf>n  rf>u   rf>df  rf>df-exp
      pi>rf  ln2>rf  euler>rf  catalan>rf  free-rfcache

      \ elementary
      -rf  |rf|
      rf+  rfn+  rfu+  rfdf+
      rf-  rfu-  urf-  rfdf-  dfrf-
      rf*  rfn*  rfu*  rf*+   rf*-
      rf/  rfn/  nrf/  rfu/   urf/  rfdf/  dfrf/
      rf-sqrt  usqrt>rf  1/rf-sqrt  rf-cbrt  rfuroot
      rf-square  rf-pow  rf^n  rf^u  u^rf  u^u>rf
      rf2^n*  rf2^u*  rf2^n/  rf2^u/
      rf-max  rf-min  rf-reldiff  rf-dim

      \ integer rounding
      rf-ceil  rf-floor  rf-trunc  rf-integer

      \ special functions
      rf-ln    rf-log2  rf-log10  rf-exp    rf-exp2  rf-exp10
      rf-cos   rf-sin   rf-tan    rf-sincos rf-sec   rf-csc   rf-cot
      rf-acos  rf-asin  rf-atan   rf-atan2
      rf-cosh  rf-sinh  rf-tanh   rf-sinhcosh
      rf-sech  rf-csch  rf-coth   rf-acosh  rf-asinh rf-atanh
      ufac>rf  rf-ln1p  rf-expm1  rf-eint   rf-li2   rf-gamma  rf-lngamma
      rf-j0  rf-j1  rf-jn  rf-y0  rf-y1  rf-yn  rf-ai
      rf-agm  rf-hypot

      \ classification
      rf-finite  rf-infinite  rf-nan  rf-regular   

      \ data manipulation
      rf-signbit  is-rf-signbit  rf-setsign  rf-copysign
      rf-nextafter  rf-nextup   rf-nextdown
      rf-checkrange rf-subnormalize

  5.9 COMPARISONS
      rfsgn
      rf>    rf=    rf<    rf>=    rf<= 
      rf0>   rf0=   rf0<   rf0>=   rf0<=
      rfn>   rfn=   rfn<   rfn>=   rfn<=
      rfu>   rfu=   rfu<   rfu>=   rfu<=
      rfdf>  rfdf=  rfdf<  rfdf>=  rfdf<=

  6 MIXED GMPFR STACKS

  6.1 BSTACK AND QSTACK
      b>q      b>qa	ba>q	  ba>qa
      b>qnum   b>qden	b>qanum	  b>qaden
      ba>qnum  ba>qden	ba>qanum  ba>qaden
      q>b      q>ba	qa>b	  qa>ba
      qnum>b   qden>b	qanum>b	  qaden>b
      qnum>ba  qden>ba	qanum>ba  qaden>ba
 
  6.2 BSTACK AND MFSTACK
      b>mf  b>mfa  ba>mf  ba>mfa
      mf>b  mf>ba  mfa>b  mfa>ba

  6.3 QSTACK AND MFSTACK
      q>mf  q>mfa  qa>mf  qa>mfa
      mf>q  mf>qa  mfa>q  mfa>qa
 
  6.4 BSTACK AND RFSTACK
      b>rf  b>rfa  ba>rf  ba>rfa
      rf>b  rf>ba  rfa>b  rfa>ba
      rf>b-exp  rf>ba-exp  rfa>b-exp  rfa>ba-exp

  6.5 QSTACK AND RFSTACK
      q>rf  q>rfa  qa>rf  qa>rfa

  6.6 MFSTACK AND RFSTACK
      mf>rf  mf>rfa  mfa>rf  mfa>rfa
      rf>mf  rf>mfa  rfa>mf  rfa>mfa

  7 RANDOM NUMBERS

  7.1 STATE
      RANDSTATE:  free-randstate  free-randstates  copy-rand

  7.2 INITIALIZATION
      randinit-default  randinit-mt  randinit-lc-2^u  randinit-lc-bits

  7.3 SEEDING
      b>seed  u>seed

  7.4 GENERATION
      rand-#bits>u    rand[0,n-1]>u
      rand-#bits>ba   rand-#bits>b   rand2-#bits>ba   rand2-#bits>b
      rand-sbits>mfa  rand-sbits>mf  rand2-limbs>mfa  rand2-limbs>mf
      rand>rfa  rand>rf  normalized-rand>rfa  normalized-rand>rf

  8 ACCUMULATOR RELEASE
      free-accs 
---


\ *** 0 NOTATION AND TERMINOLOGY
---
There are three mathematical types, or objects, in gmp,
integers, rationals (quotients), and floating-point numbers,
with C names mpz_t, mpq_t, and mpf_t.  All three are used in
mpfr, which has an additional type for reliable floating-point
numbers, mpfr_t.

Each of these corresponds to a "record" in the terminology here,
with a garbage-collected space for each that includes a stack of
record addresses.  There is a fifth data type, gmp_randstate_t,
used by both gmp and mpfr, whose references do not have a
separate stack.

gmp:    GNU Multiple Precision library.
mpfr:   GNU Multiple Precision Floating-Point Reliable library.
gmpfr:  A contraction of gmp and mpfr, used to refer to the
  Forth implementation, called forth-gmpfr.

b, q, mf, rf:  Integer big number, or just "big number",
  rational number with big-number numerator and denominator,
  multiprecision floating-point (multifloat) number,
  multiprecision reliable floating-point number, respectively.

brec, qrec, mfrec, rfrec:  A b, q, mf, or rf record, which
  manages information about the storage of the corresponding
  numerical data.

randstate:  A record that manages random number generation
  algorithms and data, corresponding to the C type
  gmp_randstate_t.

'b, 'q, 'mf, 'rf, 'rand:  The address of a brec, qrec, mfrec,
  rfrec, or randstate.  Sometimes the abbreviation 'f is used
  for 'mf or 'rf when the context is clear.

dbrec, dqrec, dmfrec, drfrec:  A brec, qrec, mfrec, or rfrec in
  bspace, qspace, mfspace, or rfspace, which means it is subject
  to garbage collection.  The "d" stands for "dynamic".  A record
  outside of all of the four spaces is called "nondynamic" or
  "external", and is not garbage collected.  Note that nondynamic
  records nevertheless manage memory that is dynamically allocated by
  gmp or mpfr.

bspace, qspace, mfspace, rfspace:  Contiguous blocks of memory
  which can be allocated by the user (words provided below).
  Each contains a header and a buffer.  The buffer contains
  dynamic records building up from lower memory and a stack of
  dynamic and nondynamic record addresses building down from
  upper memory.  The names are also used for the values BSPACE,
  QSPACE, MFSPACE and RFSPACE, which hold the addresses of
  currently active spaces.

bstack, qstack, mfstack, rfstack:  The stack in the
  corresponding space.  Although items on the stacks are
  addresses of records, their stack notation uses the style b,
  q, and f instead of 'b, 'q, and 'f.

accumulator:  A special, nondynamic brec, qrec, mfrec, or rfrec
  used for sequential calculations and/or their results.
  Currently there are BA, QA, MFA, and RFA, plus two extra
  big-number accumulators, BY, and BZ, for secondary results
  such as remainders, and one extra rf accumulator, RFX, for
  double returns such as sincos.  Accumulator arithmetic words
  have abbreviated stack effects, where what would be the
  leftmost stack input is instead in the accumulator, and what
  would be the stack output (sometimes two items) instead
  remains in the accumulator(s).

[predicate]:  Forth true/false flag corresponding to a predicate. 
  E.g., [b1<b2].

fstring:  The ANS Forth pair representation ( c-addr u) of a
  string.

s or .s suffix:  Short for an fstring, especially on the data
stack.

zstring:  A C-style sequence of characters, terminated by a null
  character.

mstring or measured string:  A counted string whose count field
  is an aligned cell instead of a character.

mzstring:  A measured string followed by a null terminator, not
  included in the count.  Strings used for multiprecision input
  and output in gmpfr are stored on the Forth side as mzstrings,
  typically appearing on the data stack in their fstring
  representation, referred to as mz.s.
  
"free":  Used in the names of words that release dynamic memory
  allocated by gmp or mpfr, and managed by a brec, qrec, mfrec,
  rfrec, or randstate.  This is aside from garbage collection
  for dynamic records.  The gmp and mpfr function names use
  "clear" instead of "free".  The use of "free" aims to avoid
  the connotation of setting to zero, which is done instead by
  gmp functions and words with "init" in their names.  For mpfr
  functions, "init" can mean setting to a NaN.

t:  Ternary return for mpfr reliable float operations, which are
  correctly rounded.  A zero t means the result is exact, taking
  the inputs to be exact, and a positive or negative t means the
  rounded result is above or below the exact result.  The
  corresponding gmpfr operations store t in the variable RFRET.
---


\ *** 1 GENERAL SUPPORT

[UNDEFINED] -rot  [IF] : -rot rot rot ; [THEN]
[UNDEFINED] rdrop [IF] : rdrop postpone r> postpone drop ; immediate [THEN]
[UNDEFINED] (f:   [IF] : (f: postpone ( ; immediate [THEN]

---
The values MZBUF and /MZBUF hold the address and size of a
buffer for transient, once only, multiple precision string input
and output for gmpfr mathematical objects.  Strings in the
buffer are mzstrings, hence the MZ in the value names.  With
some restrictions, input strings in gmp may contain embedded
whitespace characters (see the manual), including end of line
characters.  In mpfr, only leading whitespace is allowed.  The
word -MZ-WS below works around that.
---

1024 constant DEFAULT-/MZBUF
create DEFAULT-MZBUF DEFAULT-/MZBUF allot
DEFAULT-MZBUF value MZBUF    DEFAULT-/MZBUF value /MZBUF

: ?MZBUF ( u --  ) /MZBUF u> ABORT" ***MZBUF overflow" ;

: [MZ"]  \ compile:  ( "<lines><">" -- )
         \ run:      ( -- s )
  ---
    A compile-time factor used by MZ", for string input across
    lines with null termination.
  ---
  state @ 0= IF -14 throw THEN
  postpone AHEAD [char] " upto-m,s null-m+,
  2>r align postpone THEN 2r> postpone 2literal ;

: MZ" \ compile:    ( "<lines><">" -- )
      \ run:	    ( -- s )
      \ interpret:  ( "<lines><">" -- s )
  STATE @ IF
    [MZ"]
  ELSE
    MZBUF /MZBUF m!" MZBUF /MZBUF null-m+ and
    0= ABORT" ***MZBUF overflow"
    MZBUF mcount
  THEN ; immediate

[UNDEFINED] ZCOUNT [IF]
: ZCOUNT  ( 'zstr -- 'zstr len )
  dup 0 2>r
  BEGIN
    dup c@ 0<>
  WHILE
    1+ r> 1+ >r
  REPEAT drop 2r> ;
[THEN]

[UNDEFINED] ZSTRLEN [IF] : ZSTRLEN  ZCOUNT nip ; [THEN]

0 [IF]  \ both versions tested
: -mz-ws  ( mz.s -- )
  ---
    Compress mz.s in place to a new mzstring by removing all
    whitespace characters.  This implementation does one
    redundant cat if there is no leading whitespace.
  ---
  2dup -mcount ( m) >r trim| r@ off
  BEGIN |trim bl separate -rot r@ <m+> 0= UNTIL r> <null-m+> ;
[ELSE]
: -mz-ws  ( mz.s -- )
  ---
    Compress mz.s in place to a new mzstring by removing all
    whitespace characters.  This implementation does no
    redundant cats.
  ---
  2dup -mcount ( m) >r
  trim| dup 0= IF 2drop r> dup off <null-m+> EXIT THEN
  over c@ white?
  IF |trim bl separate
    IF 2drop dup negate r@ +! ELSE r@ off THEN
  ELSE r@ off THEN    
  BEGIN |trim bl separate -rot r@ <m+> 0= UNTIL r> <null-m+> ;
[THEN]


MULTI-INTEGERS [IF]
\ *** 2 INTEGERS

\ *** 2.1 RECORD SPACE
---
Since bspace is reserved with ALLOCATE, it can be released with
FREE.  However, that should be done only after executing
CLEAR-BSPACE, which is what FREE-BSPACE does; see below.
---

: copy-brec ( 'b.src 'b.targ -- )
  ---
    Assume that target record memory has been reserved. 
    Initialize it and set it to the value of the source b. 
    Intended only for internal use by gcrstack.fs.
  ---
  swap mpz_init_set ;

0 value BSPACE

: make-bspace  ( /space -- )
  ---
    This word writes over the value BSPACE.  Be sure that it's
    either saved or no longer needed.
  ---
  ['] mpz_clear ['] copy-brec rot /MPZ swap
  make-gspace to bspace ;

/MPZ 2 cells + 4000 * CONSTANT DEFAULT-/BSPACE
DEFAULT-/BSPACE make-bspace

: clear-bspace ( -- )
  ---
    Free any unreleased records in bspace, zero any bvariables
    pointing to them, and clear bspace, including the bstack.
  ---
  bspace free-grecs ;

: free-bspace ( -- ) clear-bspace bspace free throw ; 

\ *** 2.2 GARBAGE COLLECTION

: bgarbage?	   ( -- flag ) bspace ggarbage? ;
: collect-bgarbage ( -- flag ) bspace collect-ggarbage ;
: bgc-off	   ( -- )      bspace ggc-off ;
: bgc-on	   ( -- )      bspace ggc-on ;
: bgc-lock@	   ( -- flag ) bspace ggc-lock@ ;
: bgc-lock!	   ( flag -- ) bspace ggc-lock! ;

\ *** 2.3 STACK

PFE-HOST [IF]
  synonym (b: (
  synonym BVARIABLE GVARIABLE
[THEN]
GFORTH-HOST [IF]
  ' (         alias (b: immediate
  ' GVARIABLE alias BVARIABLE
[THEN]
IFORTH-HOST [IF]
  : (b:  postpone ( ; immediate
  : BVARIABLE  GVARIABLE ;
[THEN]

: bunused ( -- u )       bspace gunused ;
: bdepth  ( -- u )	 bspace gdepth ;
: (>bs)   ( 'b -- b: b ) bspace (>gs) ;
: (bdrop) (b: b -- )	 bspace (gdrop) ;

: bs@	   (b: b -- b s: bra )			     bspace gs@ ;
: bs2@	   (b: b1 b2 -- b1 b2 s: 'b1 'b2 )	     bspace gs2@ ;
: bs3@	   (b: b1 b2 b3 -- b1 b2 b3 s: 'b1 'b2 'b3 ) bspace gs3@ ;
: >bs	   ( 'b -- b: b )			     bspace >gs ;
: >bs-copy ( 'b -- b: b )			     bspace >gs-copy ;

: num-xt>bs  ( <num> xt -- b: num )
  ---
    Use xt to transfer <num>, a number on the data,
    floating-point, q, or mf stack, or a bnumber string in the
    input stream, into a dbrec, and push it onto the bstack. 
    The stack effect for xt is ( <num> 'b -- ).
  ---
  >r
  bspace >/grec @ cell+ bspace ?GROOM
  bspace new-grec ( 'b) >r
  r@ mpz_init
  2r@ swap EXECUTE
  r@ bspace (gpush)
  bspace >gsp @ r> cell- !  \ link brec to bstack
  rdrop ;

\ private
: b">'b  ( "<chars><">" 'b -- )  \ conversion defends on BASE
  dup postpone mz" ( 'b s) drop base @ mpz_set_str
  IF ( 'b) dup mpz_clear
           cell- bspace >gbrkp !    \ restore bspace state 
     true ABORT" ***invalid gmp integer string" THEN
  ( 'b) drop ;

: b" ( "<chars><">" -- b: b ) ['] b">'b num-xt>bs ;

0 value LAST-BCONSTANTP
: BCONSTANT ( "name" b: b -- )
  bs@ bspace is-grec 0=
  IF bspace gs>extern >bs-copy THEN \ new ext copy
  bspace gs>extern ( ra) >r create here r> ,
  last-bconstantp , to last-bconstantp
DOES>  (b: -- b.ext )  @ (>bs) ;

: free-bconstants ( -- )
  last-bconstantp
  BEGIN  dup 
  WHILE  dup @ mpz_clear cell+ @
  REPEAT drop ; 

: b,'b	 (b: b -- s: 'b )	    bspace g,ra ;
: b@	 ( bva -- b: b )	    bspace g@ ;
: b!	 (b: b s: bva -- )	    bspace g! ;
: b,	 (b: b -- )		    here 0 , b! ;
: bdrop  (b: b -- )		    bspace gdrop ; 
: b2drop (b: b1 b2 -- )		    bspace g2drop ;
: b3drop (b: b1 b2 b3 -- )	    bspace g3drop ;
: bdup   (b: b -- b b )		    bspace gdup ;
: b2dup  (b: b1 b2 -- b1 b2 b1 b2 ) bspace g2dup ;
: bswap  (b: b1 b2 -- b2 b1 )	    bspace gswap ;
: bnip   (b: b1 b2 -- b2 )	    bspace gnip ;
: bover  (b: b1 b2 -- b1 b2 b1 )    bspace gover ;
: btuck  (b: b1 b2 -- b2 b1 b2 )    bspace gtuck ;
: brot	 (b: b1 b2 b3 -- b2 b3 b1 ) bspace grot ;
: -brot	 (b: b1 b2 b3 -- b3 b1 b2 ) bspace -grot ;
: bpick  ( u b: bu ... b0 -- bu ... b0 bu ) bspace gpick ;
: bexchange ( i j b: maxth ... minth ... -- minth ... maxth ... )
  bspace gexchange ; 

: b.  (b: b -- )
  bs@ base @
  mpz_sizeinbase [ 1 cells 2 chars + ] literal + ?MZBUF 
  MZBUF cell+ base @ bs@ mpz_get_str bdrop
  zcount dup MZBUF ! type space ;

\ *** 2.4 ACCUMULATOR ARITHMETIC
---
An accumulator can be used for a sequence of computations to
minimize bstack transfers, but especially to optimize the
initialization of brecords and the release of the memory they
manage.

Keep in mind that computation words that input only from and
output only to the bstack overwrite not only the BA accumulator,
but also BY, when a bremainder is returned, and BY and BZ, when
gcd coefficients are returned.
---

CREATE ba /MPZ allot ba mpz_init
CREATE by /MPZ allot by mpz_init
CREATE bz /MPZ allot by mpz_init

: 0>ba  ( -- )	   ba mpz_clear ba mpz_init ;
: 0>by  ( -- )	   by mpz_clear by mpz_init ;
: 0>bz  ( -- )	   bz mpz_clear bz mpz_init ;
: b>ba  (b: b -- ) ba dup mpz_clear bs@ mpz_init_set bdrop ;
: n>ba  ( n -- )   ba dup mpz_clear swap mpz_init_set_si ;
: u>ba  ( u -- )   ba dup mpz_clear swap mpz_init_set_ui ;
: df>ba ( n -- )   ba dup mpz_clear mpz_init_set_d ;

\ Assume that BA, BY, and BZ have been initialized.

: ba>b (b: -- ba ) ba >bs-copy ;
: by>b (b: -- y )  by >bs-copy ;
: bz>b (b: -- z )  bz >bs-copy ;

: ba>n	    ( -- n )	      ba mpz_get_si ;
: ba>u	    ( -- u )	      ba mpz_get_ui ;
: ba>df	    (f: -- df )	      ba mpz_get_d ;
: ba>df-exp ( 'exp -- f: df ) ba mpz_get_d_2exp ;

: ba. ( -- ) ba>b b. ;

\ elementary
: -ba	   ( -- )	   ba dup mpz_neg ;
: |ba|	   ( -- )	   ba dup mpz_abs ;
: ba+	   (b: b -- )	   ba dup bs@ mpz_add bdrop ;
: ba-	   (b: b -- )	   ba dup bs@ mpz_sub bdrop ;
: bau+	   ( u -- )	   ba dup  rot mpz_add_ui ;
: bau-	   ( u -- )	   ba dup  rot mpz_sub_ui ;
: uba-     ( u -- )	   ba dup -rot mpz_ui_sub ;
: ba*	   (b: b -- )	   ba dup bs@ mpz_mul bdrop ;
: ban*	   ( n -- )	   ba dup rot mpz_mul_si ;
: bau*	   ( u -- )	   ba dup rot mpz_mul_ui ;
: b*a+ 	   (b: b1 b2 -- )  ba bs2@ mpz_addmul b2drop ;
: bu*a+	   (b: b s: u -- ) ba bs@ rot mpz_addmul_ui bdrop ;
: b*a-	   (b: b1 b2 -- )  ba bs2@ mpz_submul b2drop ;
: bu*a-	   (b: b s: u -- ) ba bs@ rot mpz_submul_ui bdrop ;
: ba2^u* ( #bits -- )	   >r ba dup r> mpz_mul_2exp ;

\ ceiled division
: bac/      (b: den -- )   ba dup bs@ mpz_cdiv_q bdrop ;
: bacrem    (b: den -- )   ba dup bs@ mpz_cdiv_r bdrop ;
: bac/rem   (b: den -- )   ba by ba bs@ mpz_cdiv_qr bdrop ;
: bauc/     ( u -- |rem| ) ba dup rot mpz_cdiv_q_ui ;
: baucrem   ( u -- |rem| ) ba dup rot mpz_cdiv_r_ui ;
: bauc/rem  ( u -- |rem| ) >r ba by ba r> mpz_cdiv_qr_ui ;
: bacurem   ( u -- |rem| ) ba swap mpz_cdiv_ui ;  \ unfortunate name
: ba2^uc/   ( pow -- )	   ba dup rot mpz_cdiv_q_2exp ;
: ba2^ucrem ( pow -- )	   ba dup rot mpz_cdiv_r_2exp ;

\ floored division
: baf/      (b: den -- )   ba dup bs@ mpz_fdiv_q bdrop ;
: bafrem    (b: den -- )   ba dup bs@ mpz_fdiv_r bdrop ;
: baf/rem   (b: den -- )   ba by ba bs@ mpz_fdiv_qr bdrop ;
: bauf/     ( u -- |rem| ) ba dup rot mpz_fdiv_q_ui ;
: baufrem   ( u -- |rem| ) ba dup rot mpz_fdiv_r_ui ;
: bauf/rem  ( u -- |rem| ) >r ba by ba r> mpz_fdiv_qr_ui ;
: bafurem   ( u -- |rem| ) ba swap mpz_fdiv_ui ;  \ unfortunate name
: ba2^uf/   ( pow -- )	   ba dup rot mpz_fdiv_q_2exp ;
: ba2^ufrem ( pows -- )	   ba dup rot mpz_fdiv_r_2exp ;

\ symmetric division
: bas/      (b: den -- )   ba dup bs@ mpz_tdiv_q bdrop ;
: basrem    (b: den -- )   ba dup bs@ mpz_tdiv_r bdrop ;
: bas/rem   (b: den -- )   ba by ba bs@ mpz_tdiv_qr bdrop ;
: baus/     ( u -- |rem| ) ba dup rot mpz_tdiv_q_ui ;
: bausrem   ( u -- |rem| ) ba dup rot mpz_tdiv_r_ui ;
: baus/rem  ( u -- |rem| ) >r ba by ba r> mpz_tdiv_qr_ui ;
: basurem   ( u -- |rem| ) ba swap mpz_tdiv_ui ;  \ unfortunate name
: ba2^us/   ( pow -- )	   ba dup rot mpz_tdiv_q_2exp ;
: ba2^usrem ( pow -- )	   ba dup rot mpz_tdiv_r_2exp ;

\ exact division
: ba/exact  (b: den -- ) ba dup bs@ mpz_divexact bdrop ;
: bau/exact ( u -- )	 ba dup rot mpz_divexact_ui ;

\ modulo
: bamod	 (b: den -- )      ba dup bs@ mpz_mod bdrop ;
: baumod ( u -- ba.mod.u ) ba dup rot mpz_mod_ui ;

\ powers
: ba-pow-mod     (b: pow mod -- ) ba dup bs2@ mpz_powm b2drop ; 
: ba-pow-mod-sec (b: pow mod -- ) ba dup bs2@ mpz_powm_sec b2drop ;
: ba^u-mod	 ( u b: mod -- )  ba dup rot bs@ mpz_powm_ui bdrop ;
: ba^u		 ( u -- )	  ba dup rot mpz_pow_ui ;
: u^u>ba	 ( u1 u2 -- )	  0>ba ba -rot mpz_ui_pow_ui ;

\ roots
: ba-nroot     ( u -- flag ) ba dup rot mpz_root 0<> ;
: ba-nroot-rem ( u -- )	     >r ba by ba r> mpz_rootrem ;
: ba-sqrt      ( -- )	     ba dup mpz_sqrt ;
: ba-sqrt-rem  ( -- )	     ba by ba mpz_sqrtrem ; 

\ number theoretics
: ba-next-prime	( -- )	       ba dup mpz_nextprime ;
: ba-gcd	(b: b -- )     ba dup bs@ mpz_gcd bdrop ;
: bau-gcd	( u -- gcd|0 ) ba dup rot mpz_gcd_ui ; 
: ba-gcd-st	(b: b -- )     ba 0>by 0>bz by bz ba bs@ mpz_gcdext bdrop ;
: ba-gcd-s	(b: b -- )     ba 0>by by 0 ba bs@ mpz_gcdext bdrop ;
: ba-lcm	(b: b -- )     ba dup bs@ mpz_lcm bdrop ;
: bau-lcm	( u -- )       ba dup rot mpz_lcm_ui ;
: ba-mod-invert (b: mod -- s: flag ) ba dup bs@ mpz_invert bdrop 0<> ;
: ba-remove	(b: b factor -- s: #removed ) ba dup bs@ mpz_remove bdrop ;
: fac>ba	( u -- )      0>ba ba swap mpz_fac_ui ;
: ba-bin	( u -- )      ba dup rot mpz_bin_ui ;
: bin>ba	( u1 u2 -- )  0>ba ba -rot mpz_bin_uiui ;
: fib>ba	( u -- )      0>ba ba swap mpz_fib_ui ;
: fib2>ba	( u -- )      0>ba ba 0>by by rot mpz_fib2_ui ;
: luc>ba	( u -- )      0>ba ba swap mpz_lucnum_ui ;
: luc2>ba	( u -- )      0>ba ba 0>by by rot mpz_lucnum2_ui ;

\ logic and bits
: ba-and      (b: b -- )       ba dup bs@ mpz_and bdrop ;
: ba-or	      (b: b -- )       ba dup bs@ mpz_ior bdrop ;
: ba-xor      (b: b -- )       ba dup bs@ mpz_xor bdrop ;
: ba-invert   ( -- )	       ba dup mpz_com ;
: ba-bit-on   ( bit# -- )      ba swap mpz_setbit ;
: ba-bit-off  ( bit# -- )      ba swap mpz_clrbit ;
: ba-bit-flip ( bit# -- )      ba swap mpz_combit ;
: ba-bit-test ( bit# -- flag ) ba swap mpz_tstbit 0<> ;

\ *** 2.5 ARITHMETIC

\ private
: n>'b  ( n 'b -- )	  swap mpz_set_si ;
: u>'b  ( u 'b -- )	  swap mpz_set_ui ;
: df>'b (f: df s: 'b -- ) mpz_set_d ;

: 0>b  (b: -- 0 )	   ['] drop num-xt>bs ;
: n>b  ( n -- b: n )	   ['] n>'b num-xt>bs ;
: u>b  ( u -- b: u )	   ['] u>'b num-xt>bs ;
: df>b (f: df -- b: [df] ) ['] df>'b num-xt>bs ;

: b>n	   (b: b -- s: n )	   bs@ mpz_get_si  bdrop ; 
: b>u	   (b: b -- s: |b| )	   bs@ mpz_get_ui bdrop ;
: b>df	   (b: b -- f: b )	   bs@ mpz_get_d bdrop ;
: b>df-exp (b: b s: 'exp -- f: b ) bs@ mpz_get_d_2exp bdrop ;

\ elementary
: -b	(b: b -- -b )		 b>ba -ba ba>b ;
: |b|	(b: b -- |b|)		 b>ba |ba| ba>b ;
: b+	(b: b1 b2 -- b1+b2 )	 b>ba ba+ ba>b ;
: b-	(b: b1 b2 -- b1-b2 )	 b>ba ba- -ba ba>b ;
: bu+	(b: b s: u -- b: b+u )	 b>ba bau+ ba>b ;
: bu-	(b: b s: u -- b: b-u )	 b>ba bau- ba>b ;
: ub-	(b: b s: u -- b: u-b )	 b>ba uba- ba>b ;
: b*	(b: b1 b2 -- b1*b2 )	 b>ba ba*  ba>b ;
: bn*	(b: b s: n -- b: b*n )	 b>ba ban* ba>b ;
: bu*	(b: b s: u -- b: b*u )	 b>ba bau* ba>b ;
: b*+	(b: b b1 b2 -- b+b1*b2 ) 2 bpick b>ba b*a+ bdrop ba>b ;
: b*-	(b: b b1 b2 -- b-b1*b2 ) 2 bpick b>ba b*a- bdrop ba>b ;
: bu*+	(b: b b1 s:u -- b+b1*u ) bswap b>ba bu*a+ ba>b ;
: bu*-	(b: b b1 s:u -- b-b1*u ) bswap b>ba bu*a- ba>b ;
: b2^u* (b: b s: u -- b: b*2^u ) b>ba ba2^u* ba>b ;

\ ceiled division
: bc/      (b: num den -- quo )		      bswap b>ba bac/ ba>b ;
: bcrem    (b: num den -- rem )		      bswap b>ba bacrem ba>b ;
: bc/rem   (b: num den -- quo rem )	      bswap b>ba bac/rem ba>b by>b ; 
: buc/	   (b: num s: u -- |rem| b: quo )     b>ba bauc/ ba>b ;
: bucrem   (b: num s: u -- |rem| b: rem )     b>ba baucrem ba>b ;
: buc/rem  (b: num s: u -- |rem |b: quo rem ) b>ba bauc/rem ba>b by>b ;
: bcurem   (b: num s: u -- |rem| )	      b>ba bacurem ;
: b2^uc/   (b: num s: u -- b: quo )	      b>ba ba2^uc/ ba>b ;
: b2^ucrem (b: num s: u -- b: rem )	      b>ba ba2^ucrem ba>b ;

\ floored division
: bf/      (b: num den -- quo )		      bswap b>ba baf/ ba>b ;
: bfrem    (b: num den -- rem )		      bswap b>ba bafrem ba>b ;
: bf/rem   (b: num den -- quo rem )	      bswap b>ba baf/rem ba>b by>b ; 
: buf/	   (b: num s: u -- |rem| b: quo )     b>ba bauf/ ba>b ;
: bufrem   (b: num s: u -- |rem| b: rem )     b>ba baufrem ba>b ;
: buf/rem  (b: num s: u -- |rem |b: quo rem ) b>ba bauf/rem ba>b by>b ;
: bfurem   (b: num s: u -- |rem| )	      b>ba bafurem ;
: b2^uf/   (b: num s: u -- b: quo )	      b>ba ba2^uf/ ba>b ;
: b2^ufrem (b: num s: u -- b: rem )	      b>ba ba2^ufrem ba>b ;

\ symmetric division
: bs/      (b: num den -- quo )		      bswap b>ba bas/ ba>b ;
: bsrem    (b: num den -- rem )		      bswap b>ba basrem ba>b ;
: bs/rem   (b: num den -- quo rem )	      bswap b>ba bas/rem ba>b by>b ; 
: bus/	   (b: num s: u -- |rem| b: quo )     b>ba baus/ ba>b ;
: busrem   (b: num s: u -- |rem| b: rem )     b>ba bausrem ba>b ;
: bus/rem  (b: num s: u -- |rem |b: quo rem ) b>ba baus/rem ba>b by>b ;
: bsurem   (b: num s: u -- |rem| )	      b>ba basurem ;
: b2^us/   (b: num s: u -- b: quo )	      b>ba ba2^us/ ba>b ;
: b2^usrem (b: num s: u -- b: rem )	      b>ba ba2^usrem ba>b ;

\ exact division
: b/exact  (b: num den -- num/den )	  bswap b>ba ba/exact ba>b ;
: bu/exact (b: num s: den -- b: num/den ) b>ba bau/exact ba>b ;

\ modulo
: bmod  (b: num den -- mod )	bswap b>ba bamod ba>b ;
: bumod (b: num s: den -- mod )	b>ba baumod ba>b ;

\ divisibility and congruence
: b-divisible	 (b: num den -- s: flag )   bs2@ mpz_divisible_p 0<> b2drop ;
: bu-divisible	 (b: num s: den -- flag )   bs@ swap mpz_divisible_ui_p 0<> bdrop ;
: b2^u-divisible (b: num s: pow -- flag )   bs@ swap mpz_divisible_2exp_p 0<> bdrop ;
: b-congruent	 (b: num c den -- s: flag ) bs3@ mpz_congruent_p 0<> b3drop ;
: bu-congruent	 (b: num s: c den -- flag ) bs@ -rot mpz_congruent_ui_p 0<> bdrop ;
: b2^u-congruent (b: num c s: pow -- flag ) bs2@ rot mpz_congruent_2exp_p 0<> b2drop ;

\ powers
: b-pow-mod	(b: b pow mod -- b' )	  2 bpick b>ba ba-pow-mod bdrop ba>b ;
: b-pow-mod-sec (b: b pow mod -- b' )	  2 bpick b>ba ba-pow-mod-sec bdrop ba>b ;
: b^u-mod	(b: b s: u b: mod -- b' ) bswap b>ba ba^u-mod ba>b ;
: b^u		(b: b s: u -- b: b^u )	  b>ba ba^u ba>b ;
: u^u>b		( u1 u2 -- b: u1^u2 )	  u^u>ba ba>b ;

\ roots
: b-nroot	(b: b s: u -- flag b: root ) b>ba ba-nroot ba>b ;
: b-nroot-rem	(b: b s: u -- b: root rem )  b>ba ba-nroot-rem ba>b by>b ;
: b-sqrt	(b: b -- root )		     b>ba ba-sqrt ba>b ;
: b-sqrt-rem	(b: b -- root rem )	     b>ba ba-sqrt-rem ba>b by>b ; 
: b-perfect-pow	(b: b -- flag )		     bs@ mpz_perfect_power_p 0<> bdrop ;
: b-perfect-squ	(b: b -- flag )		     bs@ mpz_perfect_square_p 0<> bdrop ;

\ number theoretics
: b-prob-prime (b: b s: #reps -- 2|1|0 )    bs@ swap mpz_probab_prime_p bdrop ;
: b-next-prime (b: b -- prime )		    b>ba ba-next-prime ba>b ;
: b-gcd	       (b: b1 b2 -- gcd )	    b>ba ba-gcd ba>b ;
: bu-gcd>u     (b: b s: u -- gcd|0 )	    0 bs@ rot mpz_gcd_ui bdrop ;
: bu-gcd  (b: b s: u -- b: gcd|b s: gcd|0 ) b>ba bau-gcd ba>b ;
: b-gcd-st     (b: b1 b2 -- s t gcd )	    bswap b>ba ba-gcd-st by>b bz>b ba>b ; 
: b-gcd-s      (b: b1 b2 -- s gcd )	    bswap b>ba ba-gcd-s by>b ba>b ; 
: b-lcm	       (b: b1 b2 -- lcm )	    b>ba ba-lcm ba>b ;
: bu-lcm       (b: b s: u -- b: lcm )	    b>ba bau-lcm ba>b ;
: b-mod-invert (b: b mod -- inv|0 s: flag ) bswap b>ba ba-mod-invert ba>b ;
: b-jacobi     (b: a b -- s: jac )	    bs2@ mpz_jacobi b2drop ;
: b-legendre   (b: a b -- s: leg )	    bs2@ mpz_legendre b2drop ;
: b-kronecker  (b: a b -- s: kro )	    bs2@ mpz_kronecker b2drop ;
: bn-kronecker (b: a s: n -- kro )	    bs@ swap mpz_kronecker_si bdrop ;
: bu-kronecker (b: a s: u -- kro )	    bs@ swap mpz_kronecker_ui bdrop ;
: nb-kronecker (b: b s: n -- kro )	    bs@ mpz_si_kronecker bdrop ;
: ub-kronecker (b: b s: u -- kro )	    bs@ mpz_ui_kronecker bdrop ;
: b-remove  (b: b factor -- b' s: #removed ) bswap b>ba ba-remove ba>b ;
: fac>b	       ( u -- b: u! )		    fac>ba ba>b ; 
: bu-bin       (b: b s: u -- b: bin )	    b>ba ba-bin ba>b ;
: bin>b	       ( u1 u2 -- b: bin )	    bin>ba ba>b ;
: fib>b	       ( u -- b: fib )	            fib>ba ba>b ;
: fib2>b       ( u -- b: fib_[u-1] fib_u )  fib2>ba by>b ba>b ;
: luc>b	       ( u -- b: luc )		    luc>ba ba>b ;
: luc2>b       ( u -- b: luc_[u-1] luc_u )  luc2>ba by>b ba>b ;

\ logic and bits
: b-and	      (b: b1 b2-- )		b>ba ba-and ba>b ;
: b-or	      (b: b1 b2 -- )		b>ba ba-or ba>b ;
: b-xor	      (b: b1 b2 -- )		b>ba ba-xor ba>b ;
: b-invert    (b: b -- ~b)		b>ba ba-invert ba>b ;
: b-pop-count (b: b -- s: #bits )       bs@ mpz_popcount bdrop ;
: b-ham-dist  (b: b1 b2 -- s: #bits )   bs2@ mpz_hamdist b2drop ; 
: b-scan0     (b: b s: start# -- bit# ) bs@ swap mpz_scan0 bdrop ;
: b-scan1     (b: b s: start# -- bit# ) bs@ swap mpz_scan1 bdrop ;
: b-bit-on    (b: b s: bit# -- b: b' )  b>ba ba-bit-on ba>b ;
: b-bit-off   (b: b s: bit# -- b: b' )  b>ba ba-bit-off ba>b ;
: b-bit-flip  (b: b s: bit# -- b: b' )  b>ba ba-bit-flip ba>b ;
: b-bit-test  (b: b s: bit# -- flag )   b>ba ba-bit-test ;

\ *** 2.6 COMPARISONS

: bsgn (b: b -- s: sgn[b] ) bs@ mpz_sgn bdrop ;

: b>  (b: b1 b2 -- s: [b1>b2] )  bs2@ mpz_cmp 0>  b2drop ;
: b<  (b: b1 b2 -- s: [b1<b2] )  bs2@ mpz_cmp 0<  b2drop ;
: b=  (b: b1 b2 -- s: [b1=b2] )  bs2@ mpz_cmp 0=  b2drop ;
: b<= (b: b1 b2 -- s: [b1<=b2] ) bs2@ mpz_cmp 0<= b2drop ;
: b>= (b: b1 b2 -- s: [b1>=b2] ) bs2@ mpz_cmp 0>= b2drop ;

: b0>  (b: b -- s: [b>0] )  bs@ 0 mpz_cmp_si 0> bdrop ;
: b0=  (b: b -- s: [b=0] )  bs@ 0 mpz_cmp_si 0= bdrop ;
: b0<  (b: b -- s: [b=0] )  bs@ 0 mpz_cmp_si 0< bdrop ;
: b0>= (b: b -- s: [b>=0] ) bs@ 0 mpz_cmp_si dup 0> swap 0= or bdrop ;
: b0<= (b: b -- s: [b<=0] ) bs@ 0 mpz_cmp_si dup 0< swap 0= or bdrop ;

: bn>  (b: b s: n -- [b>n] )  bs@ swap mpz_cmp_si 0> bdrop ;
: bn=  (b: b s: n -- [b=n] )  bs@ swap mpz_cmp_si 0= bdrop ;
: bn<  (b: b s: n -- [b<n] )  bs@ swap mpz_cmp_si 0< bdrop ;
: bn>= (b: b s: n -- [b>=n] ) bs@ swap mpz_cmp_si dup 0> swap 0= or bdrop ;
: bn<= (b: b s: n -- [b<=n] ) bs@ swap mpz_cmp_si dup 0< swap 0= or bdrop ;

: bu>  (b: b s: u -- [b>u] )  bs@ swap mpz_cmp_ui 0> bdrop ;
: bu=  (b: b s: u -- [b=u] )  bs@ swap mpz_cmp_ui 0= bdrop ;
: bu<  (b: b s: u -- [b<u] )  bs@ swap mpz_cmp_ui 0< bdrop ;
: bu>= (b: b s: u -- [b>=u] ) bs@ swap mpz_cmp_ui dup 0> swap 0= or bdrop ;
: bu<= (b: b s: u -- [b<=u] ) bs@ swap mpz_cmp_ui dup 0< swap 0= or bdrop ;

: bdf>  (b: b f: df -- s: [b>df] )  bs@ mpz_cmp_d 0> bdrop ;
: bdf=  (b: b f: df -- s: [b=df] )  bs@ mpz_cmp_d 0= bdrop ;
: bdf<  (b: b f: df -- s: [b<df] )  bs@ mpz_cmp_d 0< bdrop ;
: bdf>= (b: b f: df -- s: [b>=df] ) bs@ mpz_cmp_d dup 0> swap 0= or bdrop ;
: bdf<= (b: b f: df -- s: [b<=df] ) bs@ mpz_cmp_d dup 0< swap 0= or bdrop ;

: b|>|  (b: b1 b2 -- s: [|b1|>|b2|] )  bs2@ mpz_cmpabs 0>  b2drop ;
: b|=|  (b: b1 b2 -- s: [|b1|=|b2|] )  bs2@ mpz_cmpabs 0=  b2drop ;
: b|<|  (b: b1 b2 -- s: [|b1|<|b2|] )  bs2@ mpz_cmpabs 0<  b2drop ;
: b|>=| (b: b1 b2 -- s: [|b1|>=|b2|] ) bs2@ mpz_cmpabs 0>= b2drop ;
: b|<=| (b: b1 b2 -- s: [|b1|<=|b2|] ) bs2@ mpz_cmpabs 0<= b2drop ;

: bdf|>|  (b: b f: df -- s: [|b|>|df|] )  bs@ mpz_cmpabs_d 0> bdrop ;
: bdf|=|  (b: b f: df -- s: [|b|=|df|] )  bs@ mpz_cmpabs_d 0= bdrop ;
: bdf|<|  (b: b f: df -- s: [|b|<|df|] )  bs@ mpz_cmpabs_d 0< bdrop ;
: bdf|>=| (b: b f: df -- s: [|b|>=|df|] ) bs@ mpz_cmpabs_d dup 0> swap 0= or bdrop ;
: bdf|<=| (b: b f: df -- s: [b|<|=|df|] ) bs@ mpz_cmpabs_d dup 0< swap 0= or bdrop ;

: bu|>|  (b: b s: u -- [|b|>u] )  bs@ swap mpz_cmpabs_ui 0> bdrop ;
: bu|=|  (b: b s: u -- [|b|=u] )  bs@ swap mpz_cmpabs_ui 0= bdrop ;
: bu|<|  (b: b s: u -- [|b|<u] )  bs@ swap mpz_cmpabs_ui 0< bdrop ;
: bu|>=| (b: b s: u -- [|b|>=u] ) bs@ swap mpz_cmpabs_ui dup 0> swap 0= or bdrop ;
: bu|<=| (b: b s: u -- [|b|<=u] ) bs@ swap mpz_cmpabs_ui dup 0< swap 0= or bdrop ;
[THEN] \ MULTI-INTEGERS 


MULTI-RATIONALS [IF]
\ *** 3 RATIONALS

\ *** 3.1 RECORD SPACE
---
Since qspace is reserved with ALLOCATE, it can be released with
FREE.  However, that should be done only after executing
CLEAR-QSPACE, which is what FREE-QSPACE does; see below.
---

: copy-qrec ( 'q.src 'q.targ -- )
  ---
    Assume that memory for the target record has been reserved.
    Initialize and set it to the same value as the source.
    Intended only for internal use by gcrstack.fs.
  ---
  dup mpq_init swap mpq_set ;

0 value QSPACE

: make-qspace  ( /space -- )
  ---
    This word writes over the value QSPACE.  Be sure that it's
    either saved or no longer needed.
  ---
  ['] mpq_clear ['] copy-qrec rot /MPQ swap
  make-gspace to qspace ;

/MPQ 2 cells + 4000 * CONSTANT DEFAULT-/QSPACE
DEFAULT-/QSPACE make-qspace

: clear-qspace ( -- )
  ---
    Free any unreleased records in qspace, zero any qvariables
    pointing to them, and clear qspace, including the qstack.
  ---
  qspace free-grecs ;

: free-qspace ( -- ) clear-qspace qspace free throw ; 

\ *** 3.2 GARBAGE COLLECTION

: qgarbage?	   ( -- flag ) qspace ggarbage? ;
: collect-qgarbage ( -- flag ) qspace collect-ggarbage ;
: qgc-off	   ( -- )      qspace ggc-off ;
: qgc-on	   ( -- )      qspace ggc-on ;
: qgc-lock@	   ( -- flag ) qspace ggc-lock@ ;
: qgc-lock!	   ( flag -- ) qspace ggc-lock! ;

\ *** 3.3 STACK

PFE-HOST [IF]
  synonym (q: (
  synonym QVARIABLE GVARIABLE
[THEN]
GFORTH-HOST [IF]
  ' (         alias (q: immediate
  ' GVARIABLE alias QVARIABLE
[THEN]
IFORTH-HOST [IF]
  : (q:  postpone ( ; immediate
  : QVARIABLE  GVARIABLE ;
[THEN]

: qunused ( -- u )	 qspace gunused ;
: qdepth  ( -- u )	 qspace gdepth ;
: (>qs)	  ( 'q -- q: q ) qspace (>gs) ;
: (qdrop) (q: q -- )	 qspace (gdrop) ;

: qs@	   (q: q -- q s: qra )			     qspace gs@ ;
: qs2@	   (q: q1 q2 -- q1 q2 s: 'q1 'q2 )	     qspace gs2@ ;
: qs3@	   (q: q1 q2 q3 -- q1 q2 q3 s: 'q1 'q2 'q3 ) qspace gs3@ ;
: >qs	   ( 'q -- q: q )			     qspace >gs ;
: >qs-copy ( 'q -- q: q )			     qspace >gs-copy ;

: num-xt>qs  ( <num> xt -- q: num )
  ---
    Use xt to transfer <num>, a number on the data, b, or mf
    stack, or a q string in the input stream, into a dqrec, and
    push it onto the qstack.  The stack effect for xt is
    ( <num> 'q -- ).
  ---
  >r
  qspace >/grec @ cell+ qspace ?GROOM
  qspace new-grec ( 'q) >r
  r@ mpq_init
  2r@ swap EXECUTE
  r@ qspace (gpush)
  qspace >gsp @ r> cell- !  \ link qrec to qstack
  rdrop ;

\ private
: q">'q  ( "<chars><">" 'q -- )     \ conversion depends on BASE
  dup postpone mz" ( 'q s) drop base @ mpq_set_str
  IF ( 'q) dup mpq_clear
           cell- qspace >gbrkp !    \ restore qspace state 
     true ABORT" ***invalid gmp rational string" THEN
  ( 'q) drop ;

: q" ( "<chars><">" -- q: q ) ['] q">'q num-xt>qs ;

0 value LAST-QCONSTANTP
: QCONSTANT ( "name" q: q -- )
  qs@ qspace is-grec 0=
  IF qspace gs>extern >qs-copy THEN \ new ext copy
  qspace gs>extern ( ra) >r create here r> ,
  last-qconstantp , to last-qconstantp
DOES>  (q: -- q.ext )  @ (>qs) ;

: free-qconstants ( -- )
  last-qconstantp
  BEGIN  dup 
  WHILE  dup @ mpq_clear cell+ @
  REPEAT drop ; 

: q,'q	 (q: q -- s: 'q )	    qspace g,ra ;
: q@	 ( qva -- q: q )	    qspace g@ ;
: q!	 (q: q s: qva -- )	    qspace g! ;
: q,	 (q: q -- )		    here 0 , q! ;
: qdrop  (q: q -- )		    qspace gdrop ; 
: q2drop (q: q1 q2 -- )		    qspace g2drop ;
: q3drop (q: q1 q2 q3 -- )	    qspace g3drop ;
: qdup   (q: q -- q q )		    qspace gdup ;
: q2dup  (q: q1 q2 -- q1 q2 q1 q2 ) qspace g2dup ;
: qswap  (q: q1 q2 -- q2 q1 )	    qspace gswap ;
: qnip   (q: q1 q2 -- q2 )	    qspace gnip ;
: qover  (q: q1 q2 -- q1 q2 q1 )    qspace gover ;
: qtuck  (q: q1 q2 -- q2 q1 q2 )    qspace gtuck ;
: qrot	 (q: q1 q2 q3 -- q2 q3 q1 ) qspace grot ;
: -qrot	 (q: q1 q2 q3 -- q3 q1 q2 ) qspace -grot ;
: qpick  ( u q: qu ... q0 -- qu ... q0 qu ) qspace gpick ;
: qexchange ( i j q: maxth ... minth ... -- minth ... maxth ... )
  qspace gexchange ; 

: q.  (q: q -- )
  qs@ mpq_numref base @ mpz_sizeinbase
  qs@ mpq_denref base @ mpz_sizeinbase + cell+ 3 chars + ?MZBUF 
  MZBUF cell+ base @ qs@ mpq_get_str qdrop
  zcount dup MZBUF ! type space ;

\ *** 3.4 ACCUMULATOR ARITHMETIC
---
The discussion for b's in Sec. 2.4, ACCUMULATOR ARITHMETIC,
applies to q's as well, except that there is only one
accumulator, QA.
---

CREATE qa /MPQ allot qa mpq_init

: 0>qa   ( -- )	       qa mpq_clear qa mpq_init ;
: nu>qa ( num den -- ) 0>qa qa -rot mpq_set_si ;
: uu>qa ( num den -- ) 0>qa qa -rot mpq_set_ui ;
: df>qa (f: df -- )    0>qa qa mpq_set_d ;
: q>qa  (q: q -- )     0>qa qa qs@ mpq_set qdrop ;

\ Assume that QA has been initialized.
: qa>q   (q: -- qa )  qa >qs-copy ;
: qa>df   (f: -- df ) qa mpq_get_d ;
: qacanon ( -- )      qa mpq_canonicalize ;

: qa. ( -- ) qa>q q. ;

\ elementary
: -qa	 ( -- )	      qa dup mpq_neg ;
: |qa|	 ( -- )	      qa dup mpq_abs ;
: 1/qa	 ( -- )	      qa dup mpq_inv ;
: qa+	 (q: q -- )   qa dup qs@ mpq_add qdrop ;
: qa-	 (q: q -- )   qa dup qs@ mpq_sub qdrop ;
: qa*	 (q: q -- )   qa dup qs@ mpq_mul qdrop ;
: qa/	 (q: q -- )   qa dup qs@ mpq_div qdrop ;
: qa2^u* ( #bits -- ) >r qa dup r> mpq_mul_2exp ;
: qa2^u/ ( #bits -- ) >r qa dup r> mpq_div_2exp ;

\ *** 3.5 ARITHMETIC

\ private
: nu>'q ( num den 'q -- ) -rot mpq_set_si ;
: uu>'q ( num den 'q -- ) -rot mpq_set_ui ;
: df>'q (f: df s: 'q -- ) mpq_set_d ;
: q>q'  (q: q s: 'q -- )  qs@ mpq_set qdrop ;

: 0>q  (q: -- 0 1 )	     ['] drop  num-xt>qs ;
: nu>q ( num den -- q: n/d ) ['] nu>'q num-xt>qs ;
: uu>q ( num den -- q: n/d ) ['] uu>'q num-xt>qs ;
: df>q (f: df -- q: df )     ['] df>'q num-xt>qs ;
: qcanon (q: q -- q' )	     ['] q>q'  num-xt>qs ;

: q>df (q: q -- f: df ) qs@ mpq_get_d qdrop ;

\ elementary
: -q	(q: q -- -q )			 q>qa -qa qa>q ;
: |q|	(q: q -- |q|)			 q>qa |qa| qa>q ;
: 1/q	(q: -- 1/q )	      		 q>qa 1/qa qa>q ;
: q+	(q: q1 q2 -- q1+q2 )		 q>qa qa+ qa>q ;
: q-	(q: q1 q2 -- q1-q2 )		 q>qa qa- -qa qa>q ;
: q*	(q: q1 q2 -- q1*q2 )		 q>qa qa* qa>q ;
: q/	(q: q1 q2 -- q1/q2 )		 qswap q>qa qa/ qa>q ;
: q2^u* (q: q s: #bits -- q: q*2^#bits ) q>qa qa2^u* qa>q ;
: q2^u/ (q: q s: #bits -- q: q/2^#bits ) q>qa qa2^u/ qa>q ;

\ *** 3.6 COMPARISONS

: qsgn (q: q -- s: sgn[q] ) qs@ mpq_sgn qdrop ;

: q>  (q: q1 q2 -- s: [q1>q2] )  qs2@ mpq_cmp 0>  q2drop ;
: q<  (q: q1 q2 -- s: [q1<q2] )  qs2@ mpq_cmp 0<  q2drop ;
: q=  (q: q1 q2 -- s: [q1=q2] )  qs2@ mpq_cmp 0=  q2drop ;
: q<= (q: q1 q2 -- s: [q1<=q2] ) qs2@ mpq_cmp 0<= q2drop ;
: q>= (q: q1 q2 -- s: [q1>=q2] ) qs2@ mpq_cmp 0>= q2drop ;

: q0>  (q: q -- s: [q>0] )  qs@ 0 1 mpq_cmp_si 0> qdrop ;
: q0=  (q: q -- s: [q=0] )  qs@ 0 1 mpq_cmp_si 0= qdrop ;
: q0<  (q: q -- s: [q=0] )  qs@ 0 1 mpq_cmp_si 0< qdrop ;
: q0>= (q: q -- s: [q>=0] ) qs@ 0 1 mpq_cmp_si dup 0> swap 0= or qdrop ;
: q0<= (q: q -- s: [q<=0] ) qs@ 0 1 mpq_cmp_si dup 0< swap 0= or qdrop ;

: qnu>  (q: q s: sn ud -- [q>sn/ud] )  qs@ -rot mpq_cmp_si 0> qdrop ;
: qnu=  (q: q s: sn ud -- [q=sn/ud] )  qs@ -rot mpq_cmp_si 0= qdrop ;
: qnu<  (q: q s: sn ud -- [q<sn/ud] )  qs@ -rot mpq_cmp_si 0< qdrop ;
: qnu>= (q: q s: sn ud -- [q>=sn/ud] ) qs@ -rot mpq_cmp_si dup 0> swap 0= or qdrop ;
: qnu<= (q: q s: sn ud -- [q<=sn/ud] ) qs@ -rot mpq_cmp_si dup 0< swap 0= or qdrop ;

: quu>  (q: q s: un ud -- [q>un/ud] )  qs@ -rot mpq_cmp_ui 0> qdrop ;
: quu=  (q: q s: un ud -- [q=un/ud] )  qs@ -rot mpq_cmp_ui 0= qdrop ;
: quu<  (q: q s: un ud -- [q<un/ud] )  qs@ -rot mpq_cmp_ui 0< qdrop ;
: quu>= (q: q s: un ud -- [q>=un/ud] ) qs@ -rot mpq_cmp_ui dup 0> swap 0= or qdrop ;
: quu<= (q: q s: un ud -- [q<=un/ud] ) qs@ -rot mpq_cmp_ui dup 0< swap 0= or qdrop ;
[THEN]  \ MULTI-RATIONALS


MULTI-FLOATS [IF]
\ *** 4 FLOATING POINT

\ *** 4.1 PRECISION
---
The default precision must be set before any operation that
initializes an mfrec without copying to it from an already
initialized mfrec.  It gets set at the end of this section.

The precision of a dynamic or constant mf never changes after
being set initially, while that of the accumulator can be
changed by the user.  The words that initialize records with the
default precision are MF", plus any word ending in >MFA or >MF,
excepting MFA>MF and MF>MFA, which propagate the precision of
the source record.

In the current version of GMP, limbs normally use 32 or 64 bits,
and the minimum precision for 32-bit limbs corresponds to 2
limbs, or 64 bits.
---

: set-default-mfprecision ( #bits -- ) mpf_set_default_prec ;
: get-default-mfprecision ( -- #bits ) mpf_get_default_prec ;

128 constant QUAD-MFPRECISION
QUAD-MFPRECISION set-default-mfprecision 

\ private
6 value OUTPUT-MFPRECISION

: get-output-mfprecision ( -- u ) OUTPUT-MFPRECISION ;
: set-output-mfprecision ( u -- ) to OUTPUT-MFPRECISION ;
  ---
    Get or set the number of significand digits displayed by
    MFA. and MF. in the current base.
  ---

\ *** 4.2 RECORD SPACE
---
Since mfspace is reserved with ALLOCATE, it can be released with
FREE.  However, that should be done only after executing
CLEAR-MFSPACE, which is what FREE-MFSPACE does; see below.
---

/MPF 2 cells + 4000 * CONSTANT DEFAULT-/MFSPACE

: copy-mfrec ( 'f.src 'f.targ -- )
  ---
    Assume that memory for the target record has been reserved.
    Initialize it with the source precision, and set it to the
    same value.  Intended only for internal use by gcrstack.fs.
  ---
  swap 2dup mpf_get_prec mpf_init2 mpf_set ;

0 value MFSPACE

: mfa-start-precision ( -- addr ) mfspace >gextra1 ;

: mfa-init-precision ( --  ) get-default-mfprecision mfa-start-precision ! ;

: make-mfspace  ( /space -- )
  ---
    This word writes over the value MFSPACE.  Be sure that it's
    either saved or no longer needed.  It sets the initial
    accumulator precision to the current default.
  ---
  ['] mpf_clear ['] copy-mfrec rot /MPF swap
  make-gspace to mfspace
  mfa-init-precision ;

DEFAULT-/MFSPACE make-mfspace

: clear-mfspace ( -- )
  ---
    Free any unreleased records in mfspace, zero any mfvariables
    pointing to them, and clear mfspace, including the mfstack.
  ---
  mfspace free-grecs ;

: free-mfspace ( -- ) clear-mfspace mfspace free throw ; 


\ *** 4.3 GARBAGE COLLECTION

: mfgarbage?	    ( -- flag ) mfspace ggarbage? ;
: collect-mfgarbage ( -- flag ) mfspace collect-ggarbage ;
: mfgc-off	    ( -- )      mfspace ggc-off ;
: mfgc-on	    ( -- )      mfspace ggc-on ;
: mfgc-lock@	    ( -- flag ) mfspace ggc-lock@ ;
: mfgc-lock!	    ( flag -- ) mfspace ggc-lock! ;

\ *** 4.4 STACK

PFE-HOST [IF]
  synonym (mf: (
  synonym MFVARIABLE GVARIABLE
[THEN]
GFORTH-HOST [IF]
  ' (         alias (mf: immediate
  ' GVARIABLE alias MFVARIABLE
[THEN]
IFORTH-HOST [IF]
  : (mf:  postpone ( ; immediate
  : MFVARIABLE  GVARIABLE ;
[THEN]

: mfunused ( -- u )	   mfspace gunused ;
: mfdepth  ( -- u )	   mfspace gdepth ;
: (>mfs)   ( 'f -- mf: f ) mfspace (>gs) ;
: (mfdrop) (mf: f -- )	   mfspace (gdrop) ;

: mfs@	    (mf: f -- f s: fra )		       mfspace gs@ ;
: mfs2@	    (mf: f1 f2 -- f1 f2 s: 'f1 'f2 )	       mfspace gs2@ ;
: mfs3@	    (mf: f1 f2 f3 -- f1 f2 f3 s: 'f1 'f2 'f3 ) mfspace gs3@ ;
: >mfs	    ( 'f -- mf: f )			       mfspace >gs ;
: >mfs-copy ( 'f -- mf: f )			       mfspace >gs-copy ;

: num-xt>mfs  ( <num> xt -- mf: f )
  ---
    Use xt to transfer <num>, a number on the data, floating-
    point, b, or q stack, or an mf number string in the input
    stream, into a dmfrec using the current default precision,
    and push it onto the mfstack.  The stack effect for xt is
    ( <num> 'f -- ).
  ---
  >r
  mfspace >/grec @ cell+ mfspace ?GROOM
  mfspace new-grec ( 'f) >r
  r@ mpf_init
  2r@ swap EXECUTE
  r@ mfspace (gpush)
  mfspace >gsp @ r> cell- !  \ link mfrec to mfstack
  rdrop ;

\ private
: mf">'f  ( "<chars><">" 'f -- )   \ conversion defends on BASE
  dup postpone mz" ( 'f mz.s) drop base @ mpf_set_str
  IF ( 'f) dup mpf_clear
           cell- mfspace >gbrkp !  \ restore mfspace state 
     true ABORT" ***invalid gmp float string" THEN
  ( 'f) drop ;

: mf" ( "<chars><">" -- mf: f ) ['] mf">'f num-xt>mfs ;

0 value LAST-MFCONSTANTP
: MFCONSTANT  ( "name" mf: f -- )
  mfs@ mfspace is-grec 0=
  IF mfspace gs>extern >mfs-copy THEN \ new ext copy
  mfspace gs>extern ( ra) >r create here r> ,
  last-mfconstantp , to last-mfconstantp
DOES>  (mf: -- f.ext )  @ (>mfs) ;

: free-mfconstants  ( -- )
  last-mfconstantp
  BEGIN  dup 
  WHILE  dup @ mpf_clear cell+ @
  REPEAT drop ; 

: mf,'f	  (mf: f -- s: 'f )	      mfspace g,ra ;
: mf@	  ( mfva -- mf: f )	      mfspace g@ ;
: mf!	  (mf: f s: mfva -- )	      mfspace g! ;
: mf,	  (mf: f -- )		      here 0 , mf! ;
: mfdrop  (mf: f -- )		      mfspace gdrop ; 
: mf2drop (mf: f1 f2 -- )	      mfspace g2drop ;
: mf3drop (mf: f1 f2 f3 -- )	      mfspace g3drop ;
: mfdup   (mf: f -- f f )	      mfspace gdup ;
: mf2dup  (mf: f1 f2 -- f1 f2 f1 f2 ) mfspace g2dup ;
: mfswap  (mf: f1 f2 -- f2 f1 )	      mfspace gswap ;
: mfnip   (mf: f1 f2 -- f2 )	      mfspace gnip ;
: mfover  (mf: f1 f2 -- f1 f2 f1 )    mfspace gover ;
: mftuck  (mf: f1 f2 -- f2 f1 f2 )    mfspace gtuck ;
: mfrot	  (mf: f1 f2 f3 -- f2 f3 f1 ) mfspace grot ;
: -mfrot  (mf: f1 f2 f3 -- f3 f1 f2 ) mfspace -grot ;
: mfpick  ( u mf: fu ... f0 -- fu ... f0 fu ) mfspace gpick ;
: mfexchange ( i j mf: maxth ... minth ... -- minth ... maxth ... )
  mfspace gexchange ; 

VARIABLE MFEXP

: >MFSTR  ( u 'f -- )
  over cell+ 2 chars + ?MZBUF
  2>r MZBUF cell+ MFEXP base @ 2r> mpf_get_str zstrlen MZBUF ! ;

: MFSTR.  ( -- )
  BASE @ >r
  MZBUF mcount
  over c@ [char] - =
  IF ." -0." 1 cut-first ELSE ." 0." THEN
  type
  r@ CASE
    10 OF [char] E ENDOF
    16 OF [char] p decimal ENDOF
    [char] @ swap
  ENDCASE
  emit MFEXP @ .
  r> BASE ! ;

: (mf.) (mf: f s: u -- ) mfs@ >MFSTR MFSTR. mfdrop ;
: mf.	(mf: f -- )	 get-output-mfprecision (mf.) ;

\ *** 4.5 ACCUMULATOR ARITHMETIC
---
The discussion for b's in Sec. 2.4, ACCUMULATOR ARITHMETIC,
applies to mf's as well, except that there is only one
accumulator, MFA.
---

CREATE mfa /MPF allot mfa mpf_init

: free-mfa ( -- ) mfa mfa-start-precision @ mpf_set_prec_raw mfa mpf_clear ;  

: 0>mfa   ( -- )     free-mfa mfa mpf_init mfa-init-precision ;
: n>mfa  ( n -- )    free-mfa mfa swap mpf_init_set_si mfa-init-precision ;
: u>mfa  ( u -- )    free-mfa mfa swap mpf_init_set_ui mfa-init-precision ;
: df>mfa (f: df -- ) free-mfa mfa mpf_init_set_d  mfa-init-precision ;

: mf>mfa  (mf: f -- )
  free-mfa mfa mfs@ mpf_init_set mfdrop
  mfa mpf_get_prec mfa-start-precision ! ;

\ Assume that MFA has been initialized.

: mfa>mf     (mf: -- mfa )     mfa >mfs-copy ;
: mfa>n	     ( -- n )	       mfa mpf_get_si ;
: mfa>u	     ( -- u )	       mfa mpf_get_ui ;
: mfa>df     (f: -- df )       mfa mpf_get_d ;
: mfa>df-exp ( 'exp -- f: df ) mfa mpf_get_d_2exp ;

: (mfa.) ( u -- ) mfa >MFSTR MFSTR. ;
: mfa.	 ( -- )	  get-output-mfprecision (mfa.) ;

\ Support for changing mfa precision on the fly.
: set-mfa-precision ( #bits -- ) mfa swap mpf_set_prec_raw ;
: get-mfa-precision ( -- #bits ) mfa mpf_get_prec ;

\ elementary
: -mfa	    ( -- )	  mfa dup mpf_neg ;
: |mfa|	    ( -- )	  mfa dup mpf_abs ;
: mfa+	    (mf: f -- )	  mfa dup mfs@ mpf_add mfdrop ;
: mfa-	    (mf: f -- )	  mfa dup mfs@ mpf_sub mfdrop ;
: mfau+	    ( u -- )	  mfa dup  rot mpf_add_ui ;
: mfau-	    ( u -- )	  mfa dup  rot mpf_sub_ui ;
: umfa-	    ( u -- )	  mfa dup -rot mpf_ui_sub ;
: mfa*	    (mf: f -- )	  mfa dup mfs@ mpf_mul mfdrop ;
: mfau*	    ( u -- )  	  mfa dup rot mpf_mul_ui ;
: mfa/	    (mf: f -- )	  mfa dup mfs@ mpf_div mfdrop ;
: mfau/	    ( u -- )	  mfa dup rot mpf_div_ui ;
: umfa/	    ( u -- )	  mfa swap mfa mpf_ui_div ;
: mfa-sqrt  ( -- )	  mfa dup mpf_sqrt ;
: usqrt>mfa ( u -- )	  0>mfa mfa swap mpf_sqrt_ui ;
: mfa^u	    ( u -- )	  mfa dup rot mpf_pow_ui ;
: mfa2^u*   ( u -- )	  >r mfa dup r> mpf_mul_2exp ;
: mfa2^u/   ( u -- )	  >r mfa dup r> mpf_div_2exp ;
\ : mfarel-   (mf: ref -- ) mfa dup mfs@ swap mpf_reldiff mfdrop ; 
: mfa-reldiff (mf: f -- )  mfa dup mfs@ mpf_reldiff mfdrop ; \ mfa=|mfa-f|/mfa

\ integer rounding
: mfa-ceil    ( -- )	  mfa dup mpf_ceil ;
: mfa-floor   ( -- )	  mfa dup mpf_floor ;
: mfa-trunc   ( -- )	  mfa dup mpf_trunc ;
: mfa-integer ( -- flag ) mfa mpf_integer_p 0<> ;

\ *** 4.6 ARITHMETIC

\ private
: n>'mf    ( n 'f -- )	     swap mpf_set_si ;
: u>'mf    ( u 'f -- )	     swap mpf_set_ui ;
: df>'mf   (f: df s: 'f -- ) mpf_set_d ;
: usqrt>'mf ( u 'f -- )	     swap mpf_sqrt_ui ;

: 0>mf  (mf: -- 0 )	   ['] drop  num-xt>mfs ;
: n>mf  ( n -- mf: n )	   ['] n>'mf  num-xt>mfs ;
: u>mf  ( u -- mf: u )	   ['] u>'mf  num-xt>mfs ;
: df>mf (f: df -- mf: df ) ['] df>'mf num-xt>mfs ;

: mf>n	    (mf: f -- s: n )		mfs@ mpf_get_si mfdrop ; 
: mf>u	    (mf: f -- s: |f| )		mfs@ mpf_get_ui mfdrop ;
: mf>df	    (mf: f -- f: f )		mfs@ mpf_get_d mfdrop ;
: mf>df-exp (mf: f s: 'exp -- f: frac ) mfs@ mpf_get_d_2exp mfdrop ;

\ elementary

: -mf  (mf: f -- -f )		mf>mfa -mfa  mfa>mf ;
: |mf| (mf: f -- |f|)		mf>mfa |mfa| mfa>mf ;
: mf+  (mf: f1 f2 -- f1+f2 )	mf>mfa mfa+  mfa>mf ;
: mf-  (mf: f1 f2 -- f1-f2 )	mf>mfa mfa- -mfa mfa>mf ;
: mfu+ (mf: f s: u -- mf: f+u )	mf>mfa mfau+ mfa>mf ;
: mfu- (mf: f s: u -- mf: f-u )	mf>mfa mfau- mfa>mf ;
: umf- (mf: f s: u -- mf: u-f )	mf>mfa umfa- mfa>mf ;
: mf*  (mf: f1 f2 -- f1*f2 )	mf>mfa mfa*  mfa>mf ;
: mfu* (mf: f s: u -- mf: f*u )	mf>mfa mfau* mfa>mf ;

: mf/	   (mf: f1 f2 -- f1*f2 )    mfswap mf>mfa mfa/ mfa>mf ;
: mfu/	   (mf: f s: u -- mf: f/u ) mf>mfa mfau/ mfa>mf ;
: umf/	   ( u mf: f -- u/f )	    mf>mfa umfa/ mfa>mf ;
: mf-sqrt  (mf: f -- sqrt[f] )	    mf>mfa mfa-sqrt mfa>mf ;
: usqrt>mf ( u -- mf: sqrt[u] )	    ['] usqrt>'mf num-xt>mfs ;

: mf^u	 ( u mf: f -- f^u )	    mf>mfa mfa^u mfa>mf ;
: mf2^u* (mf: f s: u -- mf: f*2^u ) mf>mfa mfa2^u* mfa>mf ;
: mf2^u/ (mf: f s: u -- mf: f*2^u ) mf>mfa mfa2^u/ mfa>mf ;
\ : mfrel- (mf: meas ref -- reldif )  mfswap mf>mfa mfarel- mfa>mf ;
: mf-reldiff  (mf: f1 f2 -- |f1-f2|/f1 )  mfswap mf>mfa mfa-reldiff mfa>mf ;

: mf-ceil    (mf: f -- ceil[f] )  mf>mfa mfa-ceil mfa>mf ;
: mf-floor   (mf: f -- floor[f] ) mf>mfa mfa-floor mfa>mf ;
: mf-trunc   (mf: f -- trunc[f] ) mf>mfa mfa-trunc mfa>mf ;
: mf-integer (mf: f -- s: flag )  mfs@ mpf_integer_p 0<> mfdrop ;

\ *** 4.7 COMPARISONS

: mfsgn (mf: f -- s: sgn[f] ) mfs@ mpf_sgn mfdrop ;

: mf>  (mf: f1 f2 -- s: [f1>f2] )  mfs2@ mpf_cmp 0>  mf2drop ;
: mf<  (mf: f1 f2 -- s: [f1<f2] )  mfs2@ mpf_cmp 0<  mf2drop ;
: mf=  (mf: f1 f2 -- s: [f1=f2] )  mfs2@ mpf_cmp 0=  mf2drop ;
: mf<= (mf: f1 f2 -- s: [f1<=f2] ) mfs2@ mpf_cmp 0<= mf2drop ;
: mf>= (mf: f1 f2 -- s: [f1>=f2] ) mfs2@ mpf_cmp 0>= mf2drop ;

: mf0>  (mf: f -- s: [f>0] )  mfs@ 0 mpf_cmp_si 0> mfdrop ;
: mf0=  (mf: f -- s: [f=0] )  mfs@ 0 mpf_cmp_si 0= mfdrop ;
: mf0<  (mf: f -- s: [f=0] )  mfs@ 0 mpf_cmp_si 0< mfdrop ;
: mf0>= (mf: f -- s: [f>=0] ) mfs@ 0 mpf_cmp_si dup 0> swap 0= or mfdrop ;
: mf0<= (mf: f -- s: [f<=0] ) mfs@ 0 mpf_cmp_si dup 0< swap 0= or mfdrop ;

: mfn>  (mf: f s: n -- [f>n] )  mfs@ swap mpf_cmp_si 0> mfdrop ;
: mfn=  (mf: f s: n -- [f=n] )  mfs@ swap mpf_cmp_si 0= mfdrop ;
: mfn<  (mf: f s: n -- [f<n] )  mfs@ swap mpf_cmp_si 0< mfdrop ;
: mfn>= (mf: f s: n -- [f>=n] ) mfs@ swap mpf_cmp_si dup 0> swap 0= or mfdrop ;
: mfn<= (mf: f s: n -- [f<=n] ) mfs@ swap mpf_cmp_si dup 0< swap 0= or mfdrop ;

: mfu>  (mf: f s: u -- [f>u] )  mfs@ swap mpf_cmp_ui 0> mfdrop ;
: mfu=  (mf: f s: u -- [f=u] )  mfs@ swap mpf_cmp_ui 0= mfdrop ;
: mfu<  (mf: f s: u -- [f<u] )  mfs@ swap mpf_cmp_ui 0< mfdrop ;
: mfu>= (mf: f s: u -- [f>=u] ) mfs@ swap mpf_cmp_ui dup 0> swap 0= or mfdrop ;
: mfu<= (mf: f s: u -- [f<=u] ) mfs@ swap mpf_cmp_ui dup 0< swap 0= or mfdrop ;

: mfdf>  (mf: f f: df -- s: [f>df] )  mfs@ mpf_cmp_d 0> mfdrop ;
: mfdf=  (mf: f f: df -- s: [f=df] )  mfs@ mpf_cmp_d 0= mfdrop ;
: mfdf<  (mf: f f: df -- s: [f<df] )  mfs@ mpf_cmp_d 0< mfdrop ;
: mfdf>= (mf: f f: df -- s: [f>=df] ) mfs@ mpf_cmp_d dup 0> swap 0= or mfdrop ;
: mfdf<= (mf: f f: df -- s: [f<=df] ) mfs@ mpf_cmp_d dup 0< swap 0= or mfdrop ;
[THEN]  \ MFSTACK-INCLUDED


RELIABLE-FLOATS [IF]
\ *** 5 RELIABLE FLOATING POINT

PFE-HOST    [IF] LOADM libmpfr-pfe [THEN]
GFORTH-HOST [IF] s" libmpfr-gforth.fs"  REQUIRED [THEN]
IFORTH-HOST [IF] s" libmpfr-iforth.frt" INCLUDED [THEN]

\ *** 5.1 PRECISION AND RANGE
---
"Precision" refers to the same thing in mpfr and IEEE 754-2008,
the number of bits in the significand of a floating-point
number.  In mpfr, each rf record has a precision.  Records can
be initialized either with an explicit or with a default
precision. The default precision is changeable.

The precision of a dynamic or constant rf never changes after
being set initially, while that of the accumulator can be
changed by the user.  The words that initialize records with the
default precision are RF", plus any word ending in >RFA or >RF,
excepting RFA>RF and RF>RFA, which instead propagate the
precision of the source record.

Binary IEEE 754-2008 floating-point formats are specified by the
precision p and radix two, maximum and minimum exponents, emax,
and emin = 1 - emax.

The correspondence between the mpfr maximum and minimum
exponents, here called rfemax and rfemin, and the IEEE
parameters, is the following:

  rfemax = emax + 1
  rfemin = emin - p + 2 = -emax - p + 3

  nonzero, finite numbers, p-bit significands:
  2^{emin-p+1}	    <= |f| <= (2 - 2^[1-p]) x 2^emax
  2^{-1} x 2^rfemin <= |f| <= (1 - 2^[-p]) x 2^rfemax
---

: set-default-rfprecision ( #bits -- ) mpfr_set_default_prec ;
: get-default-rfprecision ( -- #bits ) mpfr_get_default_prec ;

113 set-default-rfprecision  \ IEEE 754-2008 binary128

\ private
6 value OUTPUT-RFPRECISION

: get-output-rfprecision ( -- u ) OUTPUT-RFPRECISION ;
: set-output-rfprecision ( u -- ) to OUTPUT-RFPRECISION ;
  ---
    Get or set the number of significand digits displayed by
    RFA. and RF. in the current base.
  ---

: get-rfemin ( -- rfemin) mpfr_get_emin ;
: get-rfemax ( -- rfemax) mpfr_get_emax ;

: set-rfemin ( rfemin -- flag ) mpfr_set_emin 0= ; 
: set-rfemax ( rfemax -- flag ) mpfr_set_emax 0= ; 
  ---
    The flag is true iff rfemin, respectively, rfemax is
    successfully changed and in the range of the mpfr
    implementation.  It is up to the user to ensure that all
    current rf objects have exponents within the new range after
    executing the set functions, for example, by using
    RFA-CHECKRANGE.
  ---

mpfr_get_emin_min constant RFEMIN-MIN
mpfr_get_emin_max constant RFEMIN-MAX
mpfr_get_emax_min constant RFEMAX-MIN
mpfr_get_emax_max constant RFEMAX-MAX

\ IEEE 754-2008 binary128
16384 set-rfemax drop   -16493 set-rfemin drop

: set-IEEE-754 ( p emax -- )
  dup 1+ set-rfemax	  0= ABORT" ***mpfr emax out of range"
  over set-default-rfprecision
  + negate 3 + set-rfemin 0= ABORT" ***mpfr emin out of range" ;

\ *** 5.2 ROUNDING

0 value rfrnd  \ USER: Treat this as a read-only value!

: get-rfround ( -- rnd ) mpfr_get_default_rounding_mode ;
: set-rfround ( rnd -- ) dup to rfrnd mpfr_set_default_rounding_mode ;

\ IEEE 754-2008 rounding modes
MPFR_RNDN constant RF-NEAREST	\ to nearest, even ties  
MPFR_RNDU constant RF-UP	\ toward +Inf
MPFR_RNDD constant RF-DOWN	\ toward -Inf
MPFR_RNDZ constant RF-TOWARD	\ toward zero
MPFR_RNDA constant RF-AWAY	\ away from zero
---
Also in mpfr.h:
MPFR_RNDF	faithful, not implemented, place holder only
MPFR_RNDNA	to nearest, ties away from zero
---

RF-NEAREST set-rfround  \ should be superfluous

\ The ternary int returned by the last mpfr arithmetic
\ operation:
variable rfret   rfret off

: rfret! ( -- ) postpone rfret postpone ! ; immediate

variable rfsign   rfsign off  \ used by lgamma
1 cells 8 = [IF] \ 64 bits
: rfsign@  ( -- 1 | -1 )
  ---
    Assume that rfsign holds +1 or -1 as a gmpfr int less than 64 bits. 
    Assume 2's complement arithmetic, and sign extend the -1 to 64 bits.
  ---
  rfsign @ dup 1 <> or ;  \ assumes 2's complement
[ELSE] \ 32 bits
: rfsign@  ( -- 1 | -1) postpone rfsign postpone @ ; immediate
[THEN]

\ *** 5.3 EXCEPTION FLAGS

\ mpfr 3.0.1 returns distinct bits on/off
: rfunderflow ( -- 1|0 )  mpfr_underflow_p ;
: rfoverflow  ( -- 2|0 )  mpfr_overflow_p ;
: rfinvalid   ( -- 4|0 )  mpfr_nanflag_p ;
: rfinexact   ( -- 8|0 )  mpfr_inexflag_p ;
: rferange    ( -- 16|0 ) mpfr_erangeflag_p ;

: set-rfunderflow ( -- ) mpfr_set_underflow ;
: set-rfoverflow  ( -- ) mpfr_set_overflow ;
: set-rfinvalid	  ( -- ) mpfr_set_nanflag ;
: set-rfinexact	  ( -- ) mpfr_set_inexflag ;
: set-rferange	  ( -- ) mpfr_set_erangeflag ;

: clr-rfunderflow ( -- ) mpfr_clear_underflow ;
: clr-rfoverflow  ( -- ) mpfr_clear_overflow ;
: clr-rfinvalid   ( -- ) mpfr_clear_nanflag ;
: clr-rfinexact	  ( -- ) mpfr_clear_inexflag ;
: clr-rferange	  ( -- ) mpfr_clear_erangeflag ;
: clr-rfflags	  ( -- ) mpfr_clear_flags ;

\ *** 5.4 RECORD SPACE
---
Since rfspace is reserved with ALLOCATE, it can be released with
FREE.  However, that should be done only after executing
CLEAR-RFSPACE, which is what FREE-RFSPACE does; see below.
---

/MPFR 2 cells + 4000 * CONSTANT DEFAULT-/RFSPACE

: copy-rfrec ( 'f.src 'f.targ -- )
  ---
    Assume that memory for the target record has been reserved.
    Initialize the target with the source precision and the
    default rounding mode, and set it to the same value as the
    source.  Because the target and source precisions are the
    same, there is no rounding.
  ---
  swap 2dup mpfr_get_prec mpfr_init2
  RF-NEAREST mpfr_set drop ;

0 value RFSPACE

: make-rfspace  ( /space -- )
  ---
    This word writes over the value RFSPACE.  Be sure that it's
    either saved or no longer needed.  It sets the initial
    accumulator precision to the current default.
  ---
  ['] mpfr_clear ['] copy-rfrec rot /MPFR swap
  make-gspace to rfspace ;

DEFAULT-/RFSPACE make-rfspace

: clear-rfspace ( -- )
  ---
    Free any unreleased records in rfspace, zero any rfvariables
    pointing to them, and clear rfspace, including the rfstack.
  ---
  rfspace free-grecs ;

: free-rfspace ( -- ) clear-rfspace rfspace free throw ; 

\ *** 5.5 GARBAGE COLLECTION

: rfgarbage?	    ( -- flag ) rfspace ggarbage? ;
: collect-rfgarbage ( -- flag ) rfspace collect-ggarbage ;
: rfgc-off	    ( -- )      rfspace ggc-off ;
: rfgc-on	    ( -- )      rfspace ggc-on ;
: rfgc-lock@	    ( -- flag ) rfspace ggc-lock@ ;
: rfgc-lock!	    ( flag -- ) rfspace ggc-lock! ;

\ *** 5.6 STACK

PFE-HOST [IF]
  synonym (rf: (
  synonym RFVARIABLE GVARIABLE
[THEN]
GFORTH-HOST [IF]
  ' (         alias (rf: immediate
  ' GVARIABLE alias RFVARIABLE
[THEN]
IFORTH-HOST [IF]
  : (rf:  postpone ( ; immediate
  : RFVARIABLE  GVARIABLE ;
[THEN]

: rfunused ( -- u )	   rfspace gunused ;
: rfdepth  ( -- u )	   rfspace gdepth ;
: (>rfs)   ( 'f -- rf: f ) rfspace (>gs) ;
: (rfdrop) (rf: f -- )	   rfspace (gdrop) ;

: rfs@	    (rf: f -- f s: fra )		       rfspace gs@ ;
: rfs2@	    (rf: f1 f2 -- f1 f2 s: 'f1 'f2 )	       rfspace gs2@ ;
: rfs3@	    (rf: f1 f2 f3 -- f1 f2 f3 s: 'f1 'f2 'f3 ) rfspace gs3@ ;
: >rfs	    ( 'f -- rf: f )			       rfspace >gs ;
: >rfs-copy ( 'f -- rf: f )			       rfspace >gs-copy ;

: num-xt>rfs  ( <num> xt -- rf: f )
  ---
    Use xt to transfer <num>, a number on the data, floating-
    point, b, or q stack, or an rf number string in the input
    stream, into a dynamic rfrec using the current default
    precision, pushed onto the rfstack.  The stack effect for xt
    is ( <num> 'f -- ).
 ---
  >r
  rfspace >/grec @ cell+ rfspace ?GROOM
  rfspace new-grec ( 'f) >r
  r@ mpfr_init
  2r@ swap EXECUTE
  r@ rfspace (gpush)
  rfspace >gsp @ r> cell- !  \ link rfrec to rfstack
  rdrop ;

\ private
: rf">'f  ( "<chars><">" 'f -- )
  dup postpone mz" ( 'f mz.s) over -rot -mz-ws base @ rfrnd mpfr_set_str
  IF ( 'f) dup mpfr_clear
           cell- rfspace >gbrkp !  \ restore rfspace state 
     true ABORT" ***invalid mpfr float string" THEN
  ( 'f) drop ;

: rf" ( "<chars><">" -- rf: f ) ['] rf">'f num-xt>rfs ;

0 value LAST-RFCONSTANTP
: RFCONSTANT  ( "name" rf: f -- )
  rfs@ rfspace is-grec 0=
  IF rfspace gs>extern >rfs-copy THEN \ new ext copy
  rfspace gs>extern ( ra) >r create here r> ,
  last-rfconstantp , to last-rfconstantp
DOES>  (rf: -- f.ext )  @ (>rfs) ;

: free-rfconstants  ( -- )
  last-rfconstantp
  BEGIN  dup 
  WHILE  dup @ mpfr_clear cell+ @
  REPEAT drop ; 

: rf,'f	  (rf: f -- s: 'f )	      rfspace g,ra ;
: rf@	  ( rfva -- rf: f )	      rfspace g@ ;
: rf!	  (rf: f s: rfva -- )	      rfspace g! ;
: rf,	  (rf: f -- )		      here 0 , rf! ;
: rfdrop  (rf: f -- )		      rfspace gdrop ; 
: rf2drop (rf: rf1 f2 -- )	      rfspace g2drop ;
: rf3drop (rf: f1 f2 rf3 -- )	      rfspace g3drop ;
: rfdup   (rf: f -- f f )	      rfspace gdup ;
: rf2dup  (rf: f1 f2 -- f1 f2 f1 f2 ) rfspace g2dup ;
: rfswap  (rf: f1 f2 -- f2 f1 )	      rfspace gswap ;
: rfnip   (rf: f1 f2 -- f2 )	      rfspace gnip ;
: rfover  (rf: f1 f2 -- f1 f2 f1 )    rfspace gover ;
: rftuck  (rf: f1 f2 -- f2 f1 f2 )    rfspace gtuck ;
: rfrot	  (rf: f1 f2 f3 -- f2 f3 f1 ) rfspace grot ;
: -rfrot  (rf: f1 f2 f3 -- f3 f1 f2 ) rfspace -grot ;
: rfpick  ( u rf: fu ... f0 -- fu ... f0 fu ) rfspace gpick ;
: rfexchange ( i j rf: maxth ... minth ... -- minth ... maxth ... )
  rfspace gexchange ; 

VARIABLE RFEXP

: >RFSTR  ( u 'f -- )
  over cell+ 2 chars + ?MZBUF
  2>r MZBUF cell+ RFEXP base @ 2r> rfrnd
  mpfr_get_str zstrlen MZBUF ! ;

: RFSTR.  ( -- )
  BASE @ >r
  MZBUF mcount
  over c@ [char] - =
  IF ." -0." 1 cut-first ELSE ." 0." THEN
  type
  r@ CASE
    10 OF [char] E ENDOF
    16 OF [char] p decimal RFEXP @ 4 * RFEXP ! ENDOF
    [char] @ swap
  ENDCASE
  emit RFEXP @ .
  r> BASE ! ;

: (rf.) (rf: f s: u -- ) rfs@ >RFSTR RFSTR. rfdrop ;
: rf.	(rf: f -- )	 get-output-rfprecision (rf.) ;

\ *** 5.7 ACCUMULATOR ARITHMETIC
---
The discussion for b's in Sec. 2.4, ACCUMULATOR ARITHMETIC,
applies to rf's as well, except that there are only two
accumulators, RFA and RFX.

Note that mpfr_init uses the current default precision.
---

CREATE rfa /MPFR allot rfa mpfr_init
CREATE rfx /MPFR allot rfx mpfr_init  \ for sincos, etc.

: free-rfa   ( -- ) rfa mpfr_clear ;  
: reinit-rfa ( -- ) free-rfa rfa mpfr_init ;

\ for secondary returns
: free-rfx   ( -- ) rfx mpfr_clear ;  
: reinit-rfx ( -- ) free-rfx rfx mpfr_init ; 

: set-rfa-precision ( #bits -- ) rfa swap get-rfround mpfr_prec_round rfret! ;
: get-rfa-precision ( -- #bits ) rfa mpfr_get_prec ;

:  0>rfa   ( -- ) reinit-rfa rfa  0 mpfr_set_zero ;
: -0>rfa   ( -- ) reinit-rfa rfa -1 mpfr_set_zero ;
:  nan>rfa ( -- ) reinit-rfa rfa dup false rfrnd mpfr_setsign rfret! ;
: -nan>rfa ( -- ) reinit-rfa rfa dup true  rfrnd mpfr_setsign rfret! ;
:  inf>rfa ( -- ) reinit-rfa rfa  0 mpfr_set_inf ;
: -inf>rfa ( -- ) reinit-rfa rfa -1 mpfr_set_inf ;

: n>rfa  ( n -- )    free-rfa rfa swap rfrnd mpfr_init_set_si rfret! ;
: u>rfa  ( u -- )    free-rfa rfa swap rfrnd mpfr_init_set_ui rfret! ;
: df>rfa (f: df -- ) free-rfa rfa get-rfround mpfr_init_set_d rfret! ;

: pi>rfa      ( -- ) reinit-rfa rfa rfrnd mpfr_const_pi rfret! ;
: ln2>rfa     ( -- ) reinit-rfa rfa rfrnd mpfr_const_log2 rfret! ;
: euler>rfa   ( -- ) reinit-rfa rfa rfrnd mpfr_const_euler rfret! ;
: catalan>rfa ( -- ) reinit-rfa rfa rfrnd mpfr_const_catalan rfret! ;

: rf>rfa  (rf: f -- )
  ---
    Clear and initialize RFA with the precision of f.  Set the
    RFA value to that of f, without rounding.  Note that the
    rounding mode is irrelevant because the precisions are the
    same.
  ---
  rfa rfs@ mpfr_get_prec mpfr_set_prec
  rfa rfs@ rfrnd mpfr_set drop rfdrop ; 

\ Assume that RFA and RFX have been initialized.

: rfa>rf     (rf: -- a )       rfa >rfs-copy ;
: rfa>n	     ( -- n )	       rfa rfrnd mpfr_get_si ;
: rfa>u	     ( -- u )	       rfa rfrnd mpfr_get_ui ;
: rfa>df     (f: -- df )       rfa rfrnd mpfr_get_d ;
: rfa>df-exp ( 'exp -- f: df ) rfa rfrnd mpfr_get_d_2exp ;
: rfx>rf     (rf: -- x )       rfx >rfs-copy ;

: (rfa.) ( u -- ) rfa >RFSTR RFSTR. ;
: rfa.	 ( -- )	  get-output-rfprecision  (rfa.) ;

\ elementary
: -rfa	      ( -- )	    rfa dup rfrnd mpfr_neg rfret! ;
: |rfa|	      ( -- )	    rfa dup rfrnd mpfr_abs rfret! ;
: rfa+	      (rf: f -- )   rfa dup rfs@ rfrnd mpfr_add rfret! rfdrop ;
: rfan+	      ( u -- )	    rfa dup rot  rfrnd mpfr_add_si rfret! ;
: rfau+	      ( u -- )	    rfa dup rot  rfrnd mpfr_add_ui rfret! ;
: rfadf+      (f: df -- )   rfa dup rfrnd mpfr_add_d rfret! ;
: rfa-	      (rf: f -- )   rfa dup rfs@ rfrnd mpfr_sub rfret! rfdrop ;
: rfan-	      ( n -- )	    rfa dup rot  rfrnd mpfr_sub_si rfret! ;
: rfau-	      ( u -- )	    rfa dup rot  rfrnd mpfr_sub_ui rfret! ;
: nrfa-	      ( n -- )	    rfa dup -rot rfrnd mpfr_si_sub rfret! ;
: urfa-	      ( u -- )	    rfa dup -rot rfrnd mpfr_si_sub rfret! ;
: rfadf-      (f: df -- )   rfa dup rfrnd mpfr_sub_d rfret! ;
: dfrfa-      (f: df -- )   rfa dup rfrnd mpfr_d_sub rfret! ;
: rfa*	      (rf: f -- )   rfa dup rfs@ rfrnd mpfr_mul rfret! rfdrop ;
: rfan*	      ( n -- )      rfa dup rot  rfrnd mpfr_mul_si rfret! ;
: rfau*	      ( u -- )      rfa dup rot  rfrnd mpfr_mul_ui rfret! ;
: rfadf*      (f: f -- )    rfa dup rfrnd mpfr_mul_d rfret! ;
: rfa/	      (rf: f -- )   rfa dup rfs@ rfrnd mpfr_div rfret! rfdrop ;
: rfan/	      ( n -- )	    rfa dup rot  rfrnd mpfr_div_si rfret! ;
: rfau/	      ( u -- )	    rfa dup rot  rfrnd mpfr_div_ui rfret! ;
: nrfa/	      ( n -- )	    rfa swap rfa rfrnd mpfr_si_div rfret! ;
: urfa/	      ( u -- )	    rfa swap rfa rfrnd mpfr_ui_div rfret! ;
: rfadf/      (f: df -- )   rfa dup rfrnd mpfr_div_d rfret! ;
: dfrfa/      (f: df -- )   rfa dup rfrnd mpfr_d_div rfret! ;
: rfa-sqrt    ( -- )	    rfa dup rfrnd mpfr_sqrt rfret! ;
: usqrt>rfa   ( u -- )	    reinit-rfa rfa swap rfrnd mpfr_sqrt_ui rfret! ;
: rfa-cbrt    ( -- )	    rfa dup rfrnd mpfr_cbrt rfret! ;
: rfa-uroot   ( u -- )	    rfa dup rot rfrnd mpfr_rootn_ui rfret! ;
: 1/rfa-sqrt  ( -- )	    rfa dup rfrnd mpfr_rec_sqrt rfret! ;
: rfa-square  (rf: f -- )   rfa dup rfrnd mpfr_sqr rfret! ;
: rfa-pow     (rf: f -- )   rfa dup rfs@ rfrnd mpfr_pow rfret! rfdrop ; 
: rfa^n	      ( n -- )	    rfa dup rot  rfrnd mpfr_pow_si rfret! ;
: rfa^u	      ( u -- )	    rfa dup rot  rfrnd mpfr_pow_ui rfret! ;
: u^rfa	      ( u -- )	    rfa swap over rfrnd mpfr_ui_pow rfret! ;
: u^u>rfa     ( u1 u2 -- )  reinit-rfa rfa -rot rfrnd mpfr_ui_pow_ui rfret! ;
: rfa2^n*     ( n -- )	    >r rfa dup r> rfrnd mpfr_mul_2si rfret! ;
: rfa2^u*     ( u -- )	    >r rfa dup r> rfrnd mpfr_mul_2ui rfret! ;
: rfa2^n/     ( n -- )	    >r rfa dup r> rfrnd mpfr_div_2si rfret! ;
: rfa2^u/     ( u -- )	    >r rfa dup r> rfrnd mpfr_div_2ui rfret! ;
: rfa-max     (f: f - )	    rfa dup rfs@ rfrnd mpfr_max rfret! rfdrop ;
: rfa-min     (f: f - )	    rfa dup rfs@ rfrnd mpfr_min rfret! rfdrop ;
: rfa-reldiff (rf: f -- )   rfa dup rfs@ rfrnd mpfr_reldiff rfdrop ; \ rfa=|rfa-f|/rfa
: rfa-dim     (rf: f -- )   rfa dup rfs@ rfrnd mpfr_dim rfret! rfdrop ;

\ a = a*f2 + f1
: rfa*+	    (rf: f1 f2 -- ) rfa dup rfs2@ swap rfrnd mpfr_fma rfret! rf2drop ;
\ a = a*f2 - f1
: rfa*-	    (rf: f1 f2 -- ) rfa dup rfs2@ swap rfrnd mpfr_fms rfret! rf2drop ;

\ integer rounding
: rfa-ceil    ( -- )	  rfa dup mpfr_ceil rfret! ;
: rfa-floor   ( -- )	  rfa dup mpfr_floor rfret! ;
: rfa-trunc   ( -- )	  rfa dup mpfr_trunc rfret! ;
: rfa-integer ( -- flag ) rfa mpfr_integer_p 0<> ;

\ special functions
: rfa-ln       ( -- )	   rfa dup rfrnd mpfr_log rfret! ;
: rfa-log2     ( -- )	   rfa dup rfrnd mpfr_log2 rfret! ;
: rfa-log10    ( -- )	   rfa dup rfrnd mpfr_log10 rfret! ;
: rfa-exp      ( -- )	   rfa dup rfrnd mpfr_exp rfret! ;
: rfa-exp2     ( -- )	   rfa dup rfrnd mpfr_exp2 rfret! ;
: rfa-exp10    ( -- )	   rfa dup rfrnd mpfr_exp10 rfret! ;
: rfa-cos      ( -- )	   rfa dup rfrnd mpfr_cos rfret! ;
: rfa-sin      ( -- )	   rfa dup rfrnd mpfr_sin rfret! ;
: rfa-tan      ( -- )	   rfa dup rfrnd mpfr_tan rfret! ;
: rfa-sincos   ( -- )	   rfa rfx over rfrnd mpfr_sin_cos rfret! ;
: rfa-sec      ( -- )	   rfa dup rfrnd mpfr_sec rfret! ;
: rfa-csc      ( -- )	   rfa dup rfrnd mpfr_csc rfret! ;
: rfa-cot      ( -- )	   rfa dup rfrnd mpfr_cot rfret! ;
: rfa-acos     ( -- )	   rfa dup rfrnd mpfr_acos rfret! ;
: rfa-asin     ( -- )	   rfa dup rfrnd mpfr_asin rfret! ;
: rfa-atan     ( -- )	   rfa dup rfrnd mpfr_atan rfret! ;
: rfa-atan2    (rf: y -- ) rfa rfs@ over rfrnd mpfr_atan2 rfret! rfdrop ;
: rfa-cosh     ( -- )	   rfa dup rfrnd mpfr_cosh rfret! ;
: rfa-sinh     ( -- )	   rfa dup rfrnd mpfr_sinh rfret! ;
: rfa-tanh     ( -- )	   rfa dup rfrnd mpfr_tanh rfret! ;
: rfa-sinhcosh ( -- )	   rfa rfx over rfrnd mpfr_sinh_cosh rfret! ;

: rfa-sech  ( -- )   rfa dup rfrnd mpfr_sech rfret! ;
: rfa-csch  ( -- )   rfa dup rfrnd mpfr_csch rfret! ;
: rfa-coth  ( -- )   rfa dup rfrnd mpfr_coth rfret! ;
: rfa-acosh ( -- )   rfa dup rfrnd mpfr_acosh rfret! ;
: rfa-asinh ( -- )   rfa dup rfrnd mpfr_asinh rfret! ;
: rfa-atanh ( -- )   rfa dup rfrnd mpfr_atanh rfret! ;
: ufac>rfa  ( u -- ) reinit-rfa rfa swap rfrnd mpfr_fac_ui rfret! ;

: rfa-ln1p    ( -- )   rfa dup rfrnd mpfr_log1p rfret! ;
: rfa-expm1   ( -- )   rfa dup rfrnd mpfr_expm1 rfret! ;
: rfa-eint    ( -- )   rfa dup rfrnd mpfr_eint rfret! ;
: rfa-li2     ( -- )   rfa dup rfrnd mpfr_li2 rfret! ;
: rfa-gamma   ( -- )   rfa dup rfrnd mpfr_gamma rfret! ;
: rfa-lngamma ( -- )   rfa dup rfrnd mpfr_lngamma rfret! ;
: rfa-lgamma  ( -- )   rfa rfsign over rfrnd mpfr_lgamma rfret! rfsign@ ;
: rfa-digamma ( -- )   rfa dup rfrnd mpfr_digamma rfret! ;
: rfa-zeta    ( -- )   rfa dup rfrnd mpfr_zeta rfret! ;
: uzeta>rfa   ( u -- ) reinit-rfa rfa swap rfrnd mpfr_zeta_ui rfret! ;
: rfa-erf     ( -- )   rfa dup rfrnd mpfr_erf rfret! ;
: rfa-erfc    ( -- )   rfa dup rfrnd mpfr_erfc rfret! ;
: rfa-j0      ( -- )   rfa dup rfrnd mpfr_j0 rfret! ;
: rfa-j1      ( -- )   rfa dup rfrnd mpfr_j1 rfret! ;
: rfa-jn      ( n -- ) rfa swap over rfrnd mpfr_jn rfret! ;
: rfa-y0      ( -- )   rfa dup rfrnd mpfr_y0 rfret! ;
: rfa-y1      ( -- )   rfa dup rfrnd mpfr_y1 rfret! ;
: rfa-yn      ( n -- ) rfa swap over rfrnd mpfr_yn rfret! ;
: rfa-ai      ( -- )   rfa dup rfrnd mpfr_ai rfret! ;

: rfa-agm   (rf: f -- ) rfa dup rfs@ rfrnd mpfr_agm rfret! rfdrop ;
: rfa-hypot (rf: f -- ) rfa dup rfs@ rfrnd mpfr_hypot rfret! rfdrop ;

\ classification
: rfa-finite   ( -- flag ) rfa mpfr_number_p 0<> ;
: rfa-infinite ( -- flag ) rfa mpfr_inf_p 0<> ;
: rfa-nan      ( -- flag ) rfa mpfr_nan_p 0<> ;
: rfa-regular  ( -- flag ) rfa mpfr_regular_p 0<> ;

\ data manipulation
: rfa-signbit	 ( -- minus? ) rfa mpfr_signbit 0<> ;
: rfa-setsign	 ( minus? -- ) >r rfa dup r> rfrnd mpfr_setsign rfret! ;
: rfa-copysign	 (rf: f -- )   rfa dup rfs@ rfrnd mpfr_copysign rfret! rfdrop ; 
: rfa-nextafter  (rf: f -- )   rfa rfs@ mpfr_nexttoward rfdrop ;
: rfa-nextup	 ( -- )	       rfa mpfr_nextabove ;
: rfa-nextdown	 ( -- )	       rfa mpfr_nextbelow ;
: rfa-checkrange ( t -- )      rfa swap rfrnd mpfr_check_range rfret! ;

: rfa-subnormalize ( t -- ) rfa swap rfrnd mpfr_subnormalize rfret! ;

\ *** 5.8 ARITHMETIC

\ private
: 0>'rf	      ( 'rf -- )	  0 mpfr_set_zero ;
: -0>'rf      ( 'rf -- )	 -1 mpfr_set_zero ;
: -nan>'rf    ( 'rf -- )	 dup true rfrnd mpfr_setsign drop ;
: inf>'rf     ( 'rf -- )	  0 mpfr_set_inf ;
: -inf>'rf     ( 'rf -- )	 -1 mpfr_set_inf ;
: n>'rf       ( n 'rf -- )	 swap rfrnd mpfr_set_si rfret! ;
: u>'rf       ( u 'rf -- )	 swap rfrnd mpfr_set_ui rfret! ;
: df>'rf      (f: df s: 'rf -- ) rfrnd mpfr_set_d drop ;
: usqrt>'rf   ( u 'rf -- )	 swap rfrnd mpfr_sqrt_ui rfret! ;
: u^u>'rf     ( u1 u2 'rf -- )	 -rot rfrnd mpfr_ui_pow_ui rfret! ;
: ufac>'rf    ( u 'rf -- )	 swap rfrnd mpfr_fac_ui rfret! ;
: uzeta>'rf   ( u 'rf -- )	 swap rfrnd mpfr_zeta_ui rfret! ;
: pi>'rf      ( 'rf -- )	 rfrnd mpfr_const_pi rfret! ;
: ln2>'rf     ( 'rf -- )	 rfrnd mpfr_const_log2 rfret! ;
: euler>'rf   ( 'rf -- )	 rfrnd mpfr_const_euler rfret! ;
: catalan>'rf ( 'rf -- )	 rfrnd mpfr_const_catalan rfret! ;

: free-rfcache ( -- ) mpfr_free_cache ;

: 0>rf    (rf: -- 0 )	     ['] 0>'rf    num-xt>rfs ;
: -0>rf   (rf: -- 0 )	     ['] -0>'rf   num-xt>rfs ;
: nan>rf  (rf: -- nan )	     ['] drop	  num-xt>rfs ;
: -nan>rf (rf: -- -nan )     ['] -nan>'rf num-xt>rfs ;
: inf>rf  (rf: -- +inf )     ['] inf>'rf  num-xt>rfs ;
: -inf>rf (rf: -- -nan)	     ['] -inf>'rf num-xt>rfs ;
: n>rf    ( n -- rf: n )     ['] n>'rf	  num-xt>rfs ;
: u>rf    ( u -- rf: u )     ['] u>'rf    num-xt>rfs ;
: df>rf   (f: df -- rf: df ) ['] df>'rf   num-xt>rfs ;

: pi>rf	    ( n -- rf: n )  ['] pi>'rf      num-xt>rfs ;
: ln2>rf     ( n -- rf: n ) ['] ln2>'rf	    num-xt>rfs ;
: euler>rf   ( n -- rf: n ) ['] euler>'rf   num-xt>rfs ;
: catalan>rf ( n -- rf: n ) ['] catalan>'rf num-xt>rfs ;

: rf>n	    (rf: f -- s: n )		rfs@ rfrnd mpfr_get_si rfdrop ; 
: rf>u	    (rf: f -- s: |f| )		rfs@ rfrnd mpfr_get_ui rfdrop ;
: rf>df	    (rf: f -- f: f )		rfs@ rfrnd mpfr_get_d rfdrop ;
: rf>df-exp (rf: f s: 'exp -- f: frac ) rfs@ rfrnd mpfr_get_d_2exp rfdrop ;

\ elementary

: rfop:  ( "name" xt -- )	\ name does (rf: i*f s: j*x -- rf: f' )
  ---
    Define an rftstack arithmetic operation, in terms of an RFA
    operation that returns a single result.
  ---
  >r :
  postpone rf>rfa r> compile, postpone rfa>rf
  postpone ; ;

: rfswap-op:  ( "name" xt -- ) 	\ name does (rf: i*f s: j*x -- rf: f' )
  ---
    Define an rftstack arithmetic operation with an RFSWAP
    prefix, in terms of an RFA operation that returns a single
    result.
  ---
  >r :
  postpone rfswap postpone rf>rfa r> compile, postpone rfa>rf
  postpone ; ;

' -rfa	   rfop: -rf		(rf: f -- -f )
' |rfa|	   rfop: |rf|		(rf: f -- |f|)
' rfa+	   rfop: rf+		(rf: f1 f2 -- f1+f2 )
' rfan+	   rfop: rfn+		(rf: f s: n -- rf: f+n )
' rfau+    rfop: rfu+		(rf: f s: u -- rf: f+u )
' rfadf+   rfop: rfdf+		(rf: f1 f: f2  -- rf: f1+f2 )
' rfa-	   rfswap-op: rf-	(rf: f1 f2 -- f1-f2 )
' rfan-	   rfop: rfn-		(rf: f s: n -- rf: f-n )
' rfau-	   rfop: rfu-		(rf: f s: u -- rf: f-u )
' nrfa-	   rfop: nrf-		(rf: f s: n -- rf: n-f )
' urfa-	   rfop: urf-		(rf: f s: u -- rf: u-f )
' rfadf-   rfop: rfdf-		(rf: f1 f: f2  -- rf: f1-f2 )
' dfrfa-   rfop: dfrf-		(rf: f1 f: f2  -- rf: f2-f1 )
' rfa*	   rfop: rf*		(rf: f1 f2 -- f1*f2 )
' rfan*	   rfop: rfn*		(rf: f s: n -- rf: f*n )
' rfau*	   rfop: rfu*		(rf: f s: u -- rf: f*u )
' rfa/	   rfswap-op: rf/	(rf: f1 f2 -- f1/f2 )
' rfan/	   rfop: rfn/		(rf: f s: n -- rf: f/n )
' rfau/	   rfop: rfu/		(rf: f s: u -- rf: f/u )
' nrfa/	   rfop: nrf/		( n rf: f -- n/f )
' urfa/	   rfop: urf/		( u rf: f -- u/f )
' rfadf/   rfop: rfdf/		(rf: f1 f: f2 -- rf: f1/f2 )
' dfrfa/   rfop: dfrf/		(f: f1 rf: f2 -- f1/f2 )
' rfa-sqrt rfop: rf-sqrt	(rf: f -- sqrt[f] )

: usqrt>rf ( u -- rf: sqrt[u] ) ['] usqrt>'rf num-xt>rfs ;

' 1/rfa-sqrt rfop: 1/rf-sqrt	(rf: f -- 1/sqrt[f] )
' rfa-cbrt   rfop: rf-cbrt	(rf: f -- cbrt[f] )
' rfa-uroot  rfop: rfuroot	( u rf: f -- urt[f] )
' rfa-square rfop: rf-square	(rf: f -- f*f )
' rfa-pow    rfop: rf-pow	(rf: f1 f2 -- f1^f2 )
' rfa^n	     rfop: rf^n 	( n rf: f -- f^n )
' rfa^u      rfop: rf^u 	( u rf: f -- f^u )
' u^rfa	     rfop: u^rf		( u rf: f -- u^f )

: u^u>rf ( u1 u2 -- rf: u1^u2 ) ['] u^u>'rf num-xt>rfs ;

' rfa2^n*     rfop: rf2^n*	(rf: f s: n -- rf: f*2^n )
' rfa2^u*     rfop: rf2^u*	(rf: f s: u -- rf: f*2^u )
' rfa2^n/     rfop: rf2^n/	(rf: f s: n -- rf: f*2^n )
' rfa2^u/     rfop: rf2^u/	(rf: f s: u -- rf: f*2^u )
' rfa-max     rfop: rf-max	(rf: f1 f2 -- max[f1,f2] )
' rfa-min     rfop: rf-min	(rf: f1 f2 -- min[f1,f2] )
' rfa-reldiff rfswap-op: rf-reldiff  (rf: f1 f2 -- |f1-f2|/f1 )
' rfa-dim     rfswap-op: rf-dim	(rf: f1 f2 -- dim )

' rfa*+	    rfop: rf*+		(rf: f1 f2 f3 -- f2*f3+f1 )
' rfa*-	    rfop: rf*-		(rf: f1 f2 f3 -- f2*f3-f1 )

' rfa-ceil  rfop: rf-ceil	(rf: f -- ceil[f] )
' rfa-floor rfop: rf-floor	(rf: f -- floor[f] )
' rfa-trunc rfop: rf-trunc	(rf: f -- trunc[f] )

: rf-integer (rf: f -- s: flag ) rfs@ mpfr_integer_p 0<> rfdrop ;

\ special functions

' rfa-ln    rfop: rf-ln		(rf: f -- ln[f] )
' rfa-log2  rfop: rf-log2	(rf: f -- log2[f] )
' rfa-log10 rfop: rf-log10	(rf: f -- log10[f] )
' rfa-exp   rfop: rf-exp	(rf: f -- exp[f] )
' rfa-exp2  rfop: rf-exp2	(rf: f -- exp2[f] )
' rfa-exp10 rfop: rf-exp10	(rf: f -- exp10[f] )
' rfa-cos   rfop: rf-cos	(rf: f -- cos[f] )
' rfa-sin   rfop: rf-sin	(rf: f -- sin[f] )
' rfa-tan   rfop: rf-tan	(rf: f -- tan[f] )

: rf-sincos (rf: f -- sin[f] cos[f] ) rf>rfa rfa-sincos rfa>rf rfx>rf ;

' rfa-sec   rfop: rf-sec	(rf: f -- sec[f] )
' rfa-csc   rfop: rf-csc	(rf: f -- csc[f] )
' rfa-cot   rfop: rf-cot	(rf: f -- cot[f] )
' rfa-acos  rfop: rf-acos	(rf: f -- acos[f] )
' rfa-asin  rfop: rf-asin	(rf: f -- asin[f] )
' rfa-atan  rfop: rf-atan	(rf: f -- atan[f] )
' rfa-atan2 rfop: rf-atan2	(rf: y x -- atan2[x,y] )
' rfa-cosh  rfop: rf-cosh	(rf: f -- cosh[f] )
' rfa-sinh  rfop: rf-sinh	(rf: f -- sinh[f] )
' rfa-tanh  rfop: rf-tanh	(rf: f -- tanh[f] )

: rf-sinhcosh (rf: f -- sinh[f] cosh[f] ) rf>rfa rfa-sinhcosh rfa>rf rfx>rf ;

: ufac>rf ( u -- rf: u! ) ['] ufac>'rf num-xt>rfs ;

' rfa-sech  rfop: rf-sech	(rf: f -- sech[f] )
' rfa-csch  rfop: rf-csch	(rf: f -- csch[f] )
' rfa-coth  rfop: rf-coth	(rf: f -- coth[f] )
' rfa-acosh rfop: rf-acosh	(rf: f -- acosh[f] )
' rfa-asinh rfop: rf-asinh	(rf: f -- asinh[f] )
' rfa-atanh rfop: rf-atanh	(rf: f -- atanh[f] )

' rfa-ln1p    rfop: rf-ln1p	(rf: f -- ln1p[f] )
' rfa-expm1   rfop: rf-expm1	(rf: f -- expm1[f] )
' rfa-eint    rfop: rf-eint	(rf: f -- eint[f] )
' rfa-li2     rfop: rf-li2	(rf: f -- li2[f] )
' rfa-gamma   rfop: rf-gamma	(rf: f -- gamma[f] )
' rfa-lngamma rfop: rf-lngamma	(rf: f -- lngamma[f] )
' rfa-lgamma  rfop: rf-lgamma	(rf: f -- lgamma[f] s: sign)
' rfa-digamma rfop: rf-digamma	(rf: f -- digamma[f] )
' rfa-zeta    rfop: rf-zeta	(rf: f -- zeta[f] )

: uzeta>rf ( u -- rf: zeta[u] ) ['] uzeta>'rf num-xt>rfs ;

' rfa-erf   rfop: rf-erf	(rf: f -- erf[f] )
' rfa-erfc  rfop: rf-erfc	(rf: f -- erfc[f] )
' rfa-j0    rfop: rf-j0		(rf: f -- j0[f] )
' rfa-j1    rfop: rf-j1		(rf: f -- j1[f] )
' rfa-jn    rfop: rf-jn		(rf: f s: n -- rf: jn[f] )
' rfa-y0    rfop: rf-y0		(rf: f -- y0[f] )
' rfa-y1    rfop: rf-y1		(rf: f -- y1[f] )
' rfa-yn    rfop: rf-yn		(rf: f s: n -- rf: yn[f] )
' rfa-ai    rfop: rf-ai		(rf: f -- ai[f] )
' rfa-agm   rfswap-op: rf-agm	(rf: f1 f2 -- agm[f1,f2] )
' rfa-hypot rfop: rf-hypot	(rf: f1 f2 -- hypot[f1,f2] )

\ classification
: rf-finite   (rf: f -- s: flag ) rfs@ mpfr_number_p 0<> rfdrop ;
: rf-infinite (rf: f -- s: flag ) rfs@ mpfr_inf_p 0<> rfdrop ;
: rf-nan      (rf: f -- s: flag ) rfs@ mpfr_nan_p 0<> rfdrop ;
: rf-regular  (rf: f -- s: flag ) rfs@ mpfr_regular_p 0<> rfdrop ;

\ data manipulation
: rf-signbit	(rf: f -- s: minus? )   rf>rfa rfa-signbit ;
: is-rf-signbit (rf: f -- f s: minus? ) rfs@ mpfr_signbit 0<> ;

' rfa-setsign 	   rfop:      rf-setsign      ( minus? rf: f -- f' )
' rfa-copysign	   rfswap-op: rf-copysign     (rf: f1 f2 -- f1' )
' rfa-nextafter	   rfswap-op: rf-nextafter    (rf: f1 f2 -- f1' )
' rfa-nextup	   rfop:      rf-nextup	      (rf: f -- f' )
' rfa-nextdown	   rfop:      rf-nextdown     (rf: f -- f' )
' rfa-checkrange   rfop:      rf-checkrange   ( t rf: f -- f' ) 
' rfa-subnormalize rfop:      rf-subnormalize ( t rf: f -- f' )  

\ *** 5.9 COMPARISONS

: rfsgn (rf: f -- s: sgn[f] ) rfs@ mpfr_sgn rfdrop ;

: rf>  (rf: f1 f2 -- s: [f1>f2] )  rfs2@ mpfr_greater_p 0<> rf2drop ;
: rf<  (rf: f1 f2 -- s: [f1<f2] )  rfs2@ mpfr_less_p 0<> rf2drop ;
: rf=  (rf: f1 f2 -- s: [f1=f2] )  rfs2@ mpfr_equal_p 0<> rf2drop ;
: rf<= (rf: f1 f2 -- s: [f1<=f2] ) rfs2@ mpfr_lessequal_p 0<> rf2drop ;
: rf>= (rf: f1 f2 -- s: [f1>=f2] ) rfs2@ mpfr_greaterequal_p 0<> rf2drop ;

: rf0>  (rf: f -- s: [f>0] )  rfs@ 0 mpfr_cmp_si 0> rfdrop ;
: rf0=  (rf: f -- s: [f=0] )  rfs@ mpfr_zero_p 0<> rfdrop ;
: rf0<  (rf: f -- s: [f=0] )  rfs@ 0 mpfr_cmp_si 0< rfdrop ;
: rf0>= (rf: f -- s: [f>=0] ) rfs@ 0 mpfr_cmp_si dup 0> swap 0= or rfdrop ;
: rf0<= (rf: f -- s: [f<=0] ) rfs@ 0 mpfr_cmp_si dup 0< swap 0= or rfdrop ;

: rfn>  (rf: f s: n -- [f>n] )  rfs@ swap mpfr_cmp_si 0> rfdrop ;
: rfn=  (rf: f s: n -- [f=n] )  rfs@ swap mpfr_cmp_si 0= rfdrop ;
: rfn<  (rf: f s: n -- [f<n] )  rfs@ swap mpfr_cmp_si 0< rfdrop ;
: rfn>= (rf: f s: n -- [f>=n] ) rfs@ swap mpfr_cmp_si dup 0> swap 0= or rfdrop ;
: rfn<= (rf: f s: n -- [f<=n] ) rfs@ swap mpfr_cmp_si dup 0< swap 0= or rfdrop ;

: rfu>  (rf: f s: u -- [f>u] )  rfs@ swap mpfr_cmp_ui 0> rfdrop ;
: rfu=  (rf: f s: u -- [f=u] )  rfs@ swap mpfr_cmp_ui 0= rfdrop ;
: rfu<  (rf: f s: u -- [f<u] )  rfs@ swap mpfr_cmp_ui 0< rfdrop ;
: rfu>= (rf: f s: u -- [f>=u] ) rfs@ swap mpfr_cmp_ui dup 0> swap 0= or rfdrop ;
: rfu<= (rf: f s: u -- [f<=u] ) rfs@ swap mpfr_cmp_ui dup 0< swap 0= or rfdrop ;

: rfdf>  (rf: f f: df -- s: [f>df] )  rfs@ mpfr_cmp_d 0> rfdrop ;
: rfdf=  (rf: f f: df -- s: [f=df] )  rfs@ mpfr_cmp_d 0= rfdrop ;
: rfdf<  (rf: f f: df -- s: [f<df] )  rfs@ mpfr_cmp_d 0< rfdrop ;
: rfdf>= (rf: f f: df -- s: [f>=df] ) rfs@ mpfr_cmp_d dup 0> swap 0= or rfdrop ;
: rfdf<= (rf: f f: df -- s: [f<=df] ) rfs@ mpfr_cmp_d dup 0< swap 0= or rfdrop ;
[THEN]  \ RFSTACK-INCLUDED


\ *** 6 MIXED GMPFR STACKS

MULTI-INTEGERS MULTI-RATIONALS AND [IF]
\ *** 6.1 BSTACK AND QSTACK

\ private
: b>'q  (b: b s: 'q -- ) bs@  mpq_set_z bdrop ;
: ba>'q ( 'q -- )	 ba mpq_set_z ;
: q>'b  (q: q s: 'b -- ) qs@  mpz_set_q qdrop ;
: qa>'b ( 'b -- )	 qa mpz_set_q ; 

: b>q	(b: b -- q: b )	  ['] b>'q num-xt>qs ;
: b>qa	(b: b -- )	  0>qa qa bs@  mpq_set_z bdrop ;
: ba>q	(q: -- ba )	  ['] ba>'q num-xt>qs ;
: ba>qa	(b: b -- )	  0>qa qa ba mpq_set_z ;
: q>b	(q: q -- b: [q] ) ['] q>'b num-xt>bs ;
: q>ba	(q: q -- )	  0>ba ba qs@  mpz_set_q qdrop ;
: qa>b	(b: -- [qa] )	  ['] qa>'b num-xt>bs ;
: qa>ba	(q: q -- )	  0>ba ba qa mpz_set_q ;

\ private
: b>'qn (b: num s: q' -- ) dup qs@ mpq_set qdrop bs@ mpq_set_num bdrop ;
: b>'qd (b: num s: q' -- ) dup qs@ mpq_set qdrop bs@ mpq_set_den bdrop ;
: ba>'qn (q: q s: 'q' -- ) dup qs@ mpq_set qdrop ba mpq_set_num ; 
: ba>'qd (q: q s: 'q' -- ) dup qs@ mpq_set qdrop ba mpq_set_den ; 
: qn>'b  (q: q s: 'b -- )  qs@ mpq_get_num qdrop ;
: qd>'b  (q: q s: 'b -- )  qs@ mpq_get_den qdrop ;
: qan>'b ( 'b -- )	   qa  mpq_get_num ;
: qad>'b ( 'b -- )	   qa  mpq_get_den ;

: b>qnum   (b: n' q: n/d -- n'/d ) ['] b>'qn num-xt>qs ;
: b>qden   (b: d' q: n/d -- n/d' ) ['] b>'qd num-xt>qs ;
: b>qanum  (b: n -- )		   qa bs@ mpq_set_num bdrop ;
: b>qaden  (b: d -- )		   qa bs@ mpq_set_den bdrop ;
: ba>qnum  (q: -- n'/d )	   ['] ba>'qn num-xt>qs ;
: ba>qden  (q: -- n/d' )	   ['] ba>'qd num-xt>qs ;
: ba>qanum ( -- )		   qa ba mpq_set_num ;
: ba>qaden ( -- )		   qa ba mpq_set_den ;
: qnum>b   (q: n/d -- b: n )	   ['] qn>'b num-xt>bs ;
: qden>b   (q: n/d -- b: d )	   ['] qd>'b num-xt>bs ;
: qnum>ba  (q: n/d -- )		   ba qs@ mpq_numref mpz_init_set qdrop ;
: qden>ba  (q: n/d -- )		   ba qs@ mpq_denref mpz_init_set qdrop ;
: qanum>b  ( -- num )		   ['] qan>'b num-xt>bs ;
: qaden>b  ( -- den )		   ['] qad>'b num-xt>bs ;
: qanum>ba ( -- )		   ba qa  mpq_numref mpz_init_set ;
: qaden>ba ( -- )		   ba qa  mpq_denref mpz_init_set ;
[THEN]

MULTI-INTEGERS MULTI-FLOATS AND [IF]
\ *** 6.2 BSTACK AND MFSTACK

\ private
: b>'mf	 (b: b s: 'f -- )  bs@  mpf_set_z bdrop ;
: ba>'mf	 ( 'f -- ) ba mpf_set_z ;
: mf>'b	 (mf: f s: 'b -- ) mfs@ mpz_set_f mfdrop ;
: mfa>'b ( 'b -- )	   mfa mpz_set_f ;

: b>mf	 (b: b -- mf: b )   ['] b>'mf num-xt>mfs ;
: b>mfa	 (b: b -- )	    0>mfa mfa bs@  mpf_set_z bdrop ;
: ba>mf	 (mf: -- ba )	    ['] ba>'mf num-xt>mfs ;
: ba>mfa (b: b -- )	    0>mfa mfa ba mpf_set_z ;
: mf>b	 (mf: f -- b: [f] ) ['] mf>'b num-xt>bs ;
: mf>ba	 (mf: f -- )	    0>ba ba mfs@ mpz_set_f mfdrop ;
: mfa>b	 (b: -- [mfa] )	    ['] mfa>'b num-xt>bs ;
: mfa>ba (mf: f -- )	    0>ba ba mfa mpz_set_f ;
[THEN]

MULTI-RATIONALS MULTI-FLOATS AND [IF]
\ *** 6.3 QSTACK AND MFSTACK

\ private
: q>'mf	 (q: q s: 'mf -- ) qs@  mpf_set_q qdrop ;
: qa>'mf ( 'mf -- )	   qa mpf_set_q ;
: mf>'q	 (mf: f s: 'q -- ) mfs@ mpq_set_f mfdrop ; 
: mfa>'q ( 'q -- )	   mfa mpq_set_f ;

: q>mf	 (q: q -- mf: q ) ['] q>'mf num-xt>mfs ;
: q>mfa  (q: q -- )	  0>mfa mfa qs@  mpf_set_q qdrop ;
: qa>mf	 (mf: -- qa )	  ['] qa>'mf num-xt>mfs ;
: qa>mfa (q: q -- )	  0>mfa mfa qa mpf_set_q ;
: mf>q	 (mf: f -- q: f ) ['] mf>'q num-xt>qs ;
: mf>qa	 (mf: f -- )	  0>qa qa mfs@ mpq_set_f mfdrop ;
: mfa>q  (q: -- mfa )	  ['] mfa>'q num-xt>qs ;
: mfa>qa (mf: f -- )	  0>qa qa mfa mpq_set_f ;
[THEN]

MULTI-INTEGERS RELIABLE-FLOATS AND [IF]
\ *** 6.4 BSTACK AND RFSTACK

\ private
: b>'rf	 (b: b s: 'rf -- ) bs@  rfrnd mpfr_set_z rfret! bdrop ;
: ba>'rf ( 'rf -- )	   ba   rfrnd mpfr_set_z rfret! ;
: rf>'b	 (rf: f s: 'b -- ) rfs@ rfrnd mpfr_get_z rfret! rfdrop ;
: rfa>'b ( 'b -- )	   rfa  rfrnd mpfr_get_z rfret! ;
: rf>'b-exp  (rf: f s: 'b -- ) rfs@ mpfr_get_z_2exp rfexp ! rfdrop ;
: rfa>'b-exp (rf: f s: 'b -- ) rfa  mpfr_get_z_2exp rfexp ! ;

: b>rf	 (b: b -- rf: b )   ['] b>'rf num-xt>rfs ;
: b>rfa	 (b: b -- )	    free-rfa rfa bs@ rfrnd mpfr_init_set_z rfret! bdrop ;
: ba>rf	 (rf: -- b )	    ['] ba>'rf num-xt>rfs ;
: ba>rfa ( -- )		    free-rfa rfa ba rfrnd mpfr_init_set_z rfret! ;
: rf>b	 (rf: f -- b: [f] ) ['] rf>'b num-xt>bs ;
: rf>ba	 (rf: f -- )	    0>ba ba rfs@ rfrnd mpfr_get_z rfret! rfdrop ;
: rfa>b	 (rf: f -- b: [f] ) ['] rfa>'b num-xt>bs ; 
: rfa>ba ( -- )	    	    0>ba ba rfa rfrnd mpfr_get_z rfret! ;  

: rf>b-exp   (rf: f -- b: f.sig s: exp ) ['] rf>'b-exp num-xt>bs rfexp @ ;
: rf>ba-exp  (rf: f -- s: exp )		 0>ba ba rfs@ mpfr_get_z_2exp rfdrop ;
: rfa>b-exp  ( -- b: f.sig s: exp )	 ['] rfa>'b-exp num-xt>bs rfexp @ ;
: rfa>ba-exp ( -- exp )			 0>ba ba rfa mpfr_get_z_2exp ;
[THEN]

MULTI-RATIONALS RELIABLE-FLOATS AND [IF]
\ *** 6.5 QSTACK AND RFSTACK

\ private
: q>'rf	 (q: q s: 'rf -- ) qs@ rfrnd mpfr_set_q rfret! qdrop ;
: qa>'rf ( 'rf -- )	   qa rfrnd mpfr_set_q rfret! ;

: q>rf	 (q: q -- rf: q )   ['] q>'rf num-xt>rfs ;
: q>rfa	 (q: q -- )	    free-rfa rfa qs@ rfrnd mpfr_init_set_q rfret! qdrop ;
: qa>rf	 (rf: -- q )	    ['] qa>'rf num-xt>rfs ;
: qa>rfa ( -- )		    free-rfa rfa qa rfrnd mpfr_init_set_q rfret! ;
[THEN]

MULTI-FLOATS RELIABLE-FLOATS AND [IF]
\ *** 6.6 MFSTACK AND RFSTACK

\ private
: mf>'rf  (mf: mf s: 'rf -- ) mfs@ rfrnd mpfr_set_f rfret! mfdrop ;
: mfa>'rf ( 'rf -- )	      mfa  rfrnd mpfr_set_f rfret! ;
: rf>'mf  (rf: rf s: 'mf -- ) rfs@ rfrnd mpfr_get_f rfret! rfdrop ;
: rfa>'mf ( 'mf -- )	      rfa  rfrnd mpfr_get_f rfret! ;

: mf>rf	 (mf: mf -- rf: rf ) ['] mf>'rf num-xt>rfs ;
: mf>rfa (mf: mf -- )	     free-rfa rfa mfs@ rfrnd mpfr_init_set_f rfret! mfdrop ;
: mfa>rf (rf: -- q )	     ['] mfa>'rf num-xt>rfs ;
: mfa>rfa ( -- )	     free-rfa rfa mfa rfrnd mpfr_init_set_f rfret! ;

: rf>mf	  (rf: rf -- mf: [f] ) ['] rf>'mf num-xt>mfs ;
: rf>mfa  (rf: rf -- )	       0>mfa mfa rfs@ rfrnd mpfr_get_f rfret! rfdrop ;
: rfa>mf  (rf: rf -- mf: [f] ) ['] rfa>'mf num-xt>mfs ; 
: rfa>mfa ( -- )	       0>mfa mfa rfa rfrnd mpfr_get_f rfret! ;  
[THEN]


RANDOMS [IF]
\ *** 7 RANDOM NUMBERS

\ *** 7.1 STATE
---
State words use only the data stack.
---

0 value last-randstate

: RANDSTATE:  \ compile: ( "name" -- )
              \ run:     ( -- 'rand )
  ---
    Allot an unitialized randstate and link it into the list of
    randstates.
  ---
\ create /GMP-RANDSTATE aligned allot last-randstate , ;
  create here /GMP-RANDSTATE aligned allot last-randstate ,
  ( here) to last-randstate ;

: free-randstate ( 'rand -- ) gmp_randclear ;

: free-randstates  ( -- )
  ---
    Free the randstates in the RANDVAR list.  Do not use this
    word unless all RANDVAR's are initialized.
  ---
  last-randstate
  BEGIN  dup 
  WHILE  dup gmp_randclear /GMP-RANDSTATE aligned + @
  REPEAT to last-randstate ;

: copy-rand ( 'rand.src 'rand.targ -- ) swap gmp_randinit_set ;

\ *** 7.2 INITIALIZATION
---
All words use only the data stack except for RANDINIT-LC-2^U,
which also uses the bstack.
---

\ Mersenne twister in gmp 5.0.1
: randinit-default ( 'rand -- )	gmp_randinit_default ;

\ Mersenne twister
: randinit-mt ( 'rand -- ) gmp_randinit_mt ;

MULTI-INTEGERS [IF]
\ Linear congruential:  X = (aX + c) mod 2^u
: randinit-lc-2^u  (b: a s: c u 'rand -- )
  -rot 2>r bs@ 2r> ( 'rand 'b c u) gmp_randinit_lc_2exp bdrop ;
[THEN]

\ Linear congruential:  table driven a, c, u > 2*#bits
\			u <= 128 in gmp 5.0.1
: randinit-lc-bits ( min.#bits 'rand -- flag ) swap gmp_randinit_lc_2exp_size 0<> ;

\ *** 7.3 SEEDING
---
Seeds are set from the data stack or the bstack.
---
: u>seed ( u 'rand -- )		swap gmp_randseed_ui ;
MULTI-INTEGERS [IF]
: b>seed (b: seed s: 'rand -- ) bs@ gmp_randseed bdrop ;
[THEN]

\ *** 7.4 GENERATION
---
Random numbers are generated on the data, b, mf, or rf stack, or
in the b, mf, or rf accumulator.
---

\ 0 <= rand <= 2^u - 1, 0 <= u <= bits/cell
: rand-#bits>u ( u 'rand -- rand ) swap gmp_urandomb_ui ;

MULTI-INTEGERS [IF]
\ private
: rand>'b  ( 'rand u 'b -- ) -rot mpz_urandomb ;
: rand2>'b ( 'rand u 'b -- ) -rot mpz_rrandomb ;

: rand-#bits>ba  ( u 'rand --  )	swap ba -rot mpz_urandomb ;
: rand-#bits>b	 ( u 'rand -- b: rand ) swap ['] rand>'b num-xt>bs ;
: rand2-#bits>ba ( u 'rand --  )	swap ba -rot mpz_rrandomb ;
: rand2-#bits>b  ( u 'rand -- b: rand ) swap ['] rand2>'b num-xt>bs ;
[THEN]

MULTI-FLOATS [IF]
\ private
: rand>'mf  ( 'rand sbits 'rf -- ) -rot mpf_urandomb ;
: rand2>'mf ( limbs p 'rf -- )	   -rot mpf_random2 ;

\ 0 <= rand < 1,  sbits significant bits
: rand-sbits>mfa ( sbits 'rand -- )	     swap mfa -rot mpf_urandomb ;
: rand-sbits>mf	 ( sbits 'rand -- fs: rand ) swap ['] rand>'mf num-xt>mfs ;

\ -1.0 x lbase^-p < |rand| < 1.0 x lbase^p [inequality?!]
\	sgn[rand] = sgn[max#limbs]
\       note the absence of 'rand
: rand2-limbs>mfa ( max#limbs u -- )	       mfa -rot mpf_random2 ;
: rand2-limbs>mf  ( max#limbs u -- fs: rand ) ['] rand2>'mf num-xt>mfs ;
[THEN]

RELIABLE-FLOATS [IF]
\ private
: rand>'rf  ( 'rand 'rf -- ) swap rfrnd mpfr_urandom rfret! ;
: randb>'rf ( 'rand 'rf -- ) swap mpfr_urandomb
  ABORT" ***can't normalize random rf" ;

: rand>rfa	      ( 'rand -- )	    rfa rand>'rf ;
: rand>rf	      ( 'rand -- rf: rand ) ['] rand>'rf num-xt>rfs ;
: normalized-rand>rfa ( 'rand -- )	    rfa randb>'rf ;
: normalized-rand>rf  ( 'rand -- rf: rand ) ['] randb>'rf num-xt>rfs ;
[THEN]

\ 0<= rand < n,  rand and n are u's
: rand[0,n-1]>u ( +n 'rand -- rand ) swap gmp_urandomm_ui ;
[THEN]  \ RANDOMS

\ *** 8 ACCUMULATOR RELEASE

: free-accs  ( -- )
[ MULTI-INTEGERS ]  [IF] ba mpz_clear  by mpz_clear  bz mpz_clear [THEN]
[ MULTI-RATIONALS ] [IF] qa mpq_clear  [THEN]
[ MULTI-FLOATS ]    [IF] free-mfa      [THEN]
[ RELIABLE-FLOATS ] [IF] rfa mpfr_clear rfx mpfr_clear [THEN]
;