Gnu MPC interface

[Marcel Hendrix]

Krishna Myneni <krishna.myn...@ccreweb.org> writes Re: ZGAMMA revisited

> On Feb 4, 2:11 pm, m...@iae.nl (Marcel Hendrix) wrote:
>> m...@iae.nl (Marcel Hendrix) wrote Re: ZGAMMA revisited
...
>> The accuracy has increased to 1e-41. I have added Godfrey's/Toth's
>> P-generator, so it is now, in principal (*), possible to configure for almost
>> arbitrary precision.
...

> Can you post your mpc library bindings? They would be a useful
> starting point for generating bindings for other Forth systems.

These are the iForth bindings and the high-level user interface.
The bindings are 64-bit, although that shouldn't make a difference.
An 'sxint' is used unstead of an 'int' when C returns a 32-bit number
that must be sign-extended.

The functions with '_macro' are macro's in the DLL/so/dylib. I added
some code in src/clear.c (I compiled mpc. myself) to convert them to
functions. If this is not done an extra DLL/so is necessary.

I checked today that my extensive source changes to gmp, mpfr and mpc for
Windows 64 still allow a build on Linux that passes all check(s).
The results for e.g. Z#GAMMA and Z#WOFZ are bit for bit the same on
Windows and Linux.

The mpc.frt file borrows from a very similar mpfr.frt file. Where one
sees a 'Z#' prefix in mpc, the code in mpfr is almost the same but
uses an 'F#' prefix. Both mpc and mpfr feature a circular stack (no
reset needed) to keep numbers on.

The mpc interface uses mpfr to print numbers. The mpc call to do this
is very inflexible.

A Z#VALUE is IMMEDIATE. I may need to move it into the kernel to get
rid of that. Likewise, there are no Z#LOCALs yet.

Z#.S, .ZMAT, .ZARR and Z^^ are debugging tools. They might disappear
from the interface or get integrated into the developer interface at
a later stage.

All words not marked PRIVATE are globally visible.

-marcel

-- ------------------------------------------------------------------------
(*
* LANGUAGE : Forth
* PROJECT : iForth mpc binding
* DESCRIPTION : Interface for libmpc
* CATEGORY : Multiple precision
* AUTHOR : Marcel Hendrix
* STARTED : Thursday, January 12, 2012, 00:38
* CHANGED : Friday, January 27, 2012, 20:12, mhx
*)

NEEDS -miscutil
NEEDS -mpfr

REVISION -mpc "--- MPC 0.9 Version 2.01 ---"

PRIVATES

WINDOWS? [IF] Library: libmpc-2.dll [THEN]
LINUX? [IF] Library: libmpc.so [THEN]
DARWIN? [IF] Library: libmpc.dylib [THEN]

-- libmpc 0.9 functions ( a subset )

AliasedExtern: mpc_init2 void mpc_init2 (a, n); PRIVATE
AliasedExtern: mpc_set sxint mpc_set (a, a, n); PRIVATE
AliasedExtern: mpc_set_fr_fr sxint mpc_set_fr_fr (a, re, im, rnd); PRIVATE
AliasedExtern: mpc_set_si_si sxint mpc_set_si_si (a, n1, n2, rnd); PRIVATE
AliasedExtern: mpc_set_str sxint mpc_set_str (a, a, n, n); PRIVATE
AliasedExtern: mpc_get_str char* mpc_get_str (n, n, a, n); PRIVATE
AliasedExtern: mpc_free_str void mpc_free_str (a); PRIVATE
AliasedExtern: mpc_swap void mpc_swap (a, a); PRIVATE

AliasedExtern: mpc_get_prec int mpc_get_prec (a); PRIVATE
AliasedExtern: mpc_set_prec int mpc_set_prec (a, n); PRIVATE

AliasedExtern: mpc_cmp int mpc_cmp (a, a); PRIVATE

AliasedExtern: mpc_real void mpc_real (a, a, n); PRIVATE
AliasedExtern: mpc_imag void mpc_imag (a, a, n); PRIVATE
AliasedExtern: mpc_arg void mpc_arg (a, a, n); PRIVATE

AliasedExtern: mpc_add sxint mpc_add (a, a, a, n); PRIVATE
AliasedExtern: mpc_add_si sxint mpc_add_si (a, a, n, n); PRIVATE
AliasedExtern: mpc_add_fr sxint mpc_add_fr (a, a, a#f, n); PRIVATE
AliasedExtern: mpc_sub sxint mpc_sub (a, a, a, n); PRIVATE
AliasedExtern: mpc_sub_ui sxint mpc_sub_ui (a, a, u, n); PRIVATE
AliasedExtern: mpc_sub_fr sxint mpc_sub_fr (a, a, a#f, n); PRIVATE
AliasedExtern: mpc_mul sxint mpc_mul (a, a, a, n); PRIVATE
AliasedExtern: mpc_mul_si sxint mpc_mul_si (a, a, n, n); PRIVATE
AliasedExtern: mpc_mul_fr sxint mpc_mul_fr (a, a, a#f, n); PRIVATE
AliasedExtern: mpc_div sxint mpc_div (a, a, a, n); PRIVATE
AliasedExtern: mpc_div_ui sxint mpc_div_ui (a, a, u, n); PRIVATE
AliasedExtern: mpc_ui_div sxint mpc_ui_div (a, n, a, n); PRIVATE
AliasedExtern: mpc_div_fr sxint mpc_div_fr (a, a, f#, n); PRIVATE
AliasedExtern: mpc_fr_div sxint mpc_fr_div (a, f#, a, n); PRIVATE

AliasedExtern: mpc_neg sxint mpc_neg (a, a, n); PRIVATE
AliasedExtern: mpc_conj sxint mpc_conj (a, a, n); PRIVATE
AliasedExtern: mpc_abs sxint mpc_abs (a, a, n); PRIVATE
AliasedExtern: mpc_norm sxint mpc_norm (a, a, n); PRIVATE

AliasedExtern: mpc_sqr sxint mpc_sqr (a, a, n); PRIVATE
AliasedExtern: mpc_sqrt sxint mpc_sqrt (a, a, n); PRIVATE
AliasedExtern: mpc_pow sxint mpc_pow (a, a, a, n); PRIVATE
AliasedExtern: mpc_pow_si sxint mpc_pow_si (a, a, n, n); PRIVATE
PRIVATE
AliasedExtern: mpc_log sxint mpc_log (a, a, n); PRIVATE
AliasedExtern: mpc_exp sxint mpc_exp (a, a, n); PRIVATE

AliasedExtern: mpc_sin sxint mpc_sin (a, a, n); PRIVATE
AliasedExtern: mpc_cos sxint mpc_cos (a, a, n); PRIVATE
AliasedExtern: mpc_sin_cos sxint mpc_sin_cos (a, a, a, n); PRIVATE
AliasedExtern: mpc_tan sxint mpc_tan (a, a, n); PRIVATE
AliasedExtern: mpc_cosh sxint mpc_cosh (a, a, n); PRIVATE
AliasedExtern: mpc_sinh sxint mpc_sinh (a, a, n); PRIVATE
AliasedExtern: mpc_tanh sxint mpc_tanh (a, a, n); PRIVATE
AliasedExtern: mpc_acos sxint mpc_acos (a, a, n); PRIVATE
AliasedExtern: mpc_asin sxint mpc_asin (a, a, n); PRIVATE
AliasedExtern: mpc_atan sxint mpc_atan (a, a, n); PRIVATE
AliasedExtern: mpc_acosh sxint mpc_acosh (a, a, n); PRIVATE
AliasedExtern: mpc_asinh sxint mpc_asinh (a, a, n); PRIVATE
AliasedExtern: mpc_atanh sxint mpc_atanh (a, a, n); PRIVATE

AliasedExtern: mpc_realref int mpc_realref_macro (a); PRIVATE
AliasedExtern: mpc_imagref int mpc_imagref_macro (a); PRIVATE

-- Specials
AliasedExtern: MPC_RNDNN int MPC_RNDNN_macro (void);
AliasedExtern: MPC_RNDNZ int MPC_RNDNZ_macro (void);
AliasedExtern: MPC_RNDNU int MPC_RNDNU_macro (void);
AliasedExtern: MPC_RNDND int MPC_RNDND_macro (void);

AliasedExtern: MPC_RNDZN int MPC_RNDZN_macro (void);
AliasedExtern: MPC_RNDZZ int MPC_RNDZZ_macro (void);
AliasedExtern: MPC_RNDZU int MPC_RNDZU_macro (void);
AliasedExtern: MPC_RNDZD int MPC_RNDZD_macro (void);

AliasedExtern: MPC_RNDUN int MPC_RNDUN_macro (void);
AliasedExtern: MPC_RNDUZ int MPC_RNDUZ_macro (void);
AliasedExtern: MPC_RNDUU int MPC_RNDUU_macro (void);
AliasedExtern: MPC_RNDUD int MPC_RNDUD_macro (void);

AliasedExtern: MPC_RNDDN int MPC_RNDDN_macro (void);
AliasedExtern: MPC_RNDDZ int MPC_RNDDZ_macro (void);
AliasedExtern: MPC_RNDDU int MPC_RNDDU_macro (void);
AliasedExtern: MPC_RNDDD int MPC_RNDDD_macro (void);

AliasedExtern: /MPC int sizeof_mpc (void);

MPC_RNDNN VALUE Z#RND ( round to nearest )
#10 VALUE Z#BASE
GET-F#PRECISION VALUE Z#PRECISION PRIVATE

: SET-Z#PRECISION ( #bits -- ) TO Z#PRECISION ; GET-F#PRECISION SET-Z#PRECISION
: GET-Z#PRECISION ( -- #bits ) Z#PRECISION ;

Z#PRECISION 2* =: /NUMPAD PRIVATE
CREATE inumpad PRIVATE /NUMPAD CHARS ALLOT \ good to about 500 bits

-- Initialized to 0
: Z#BLOCK ( n | name -- )
CREATE HERE OVER /MPC * ALLOT
( n here -- )
SWAP 0 ?DO DUP Z#PRECISION mpc_init2
DUP 0 0 Z#RND mpc_set_si_si DROP
/MPC +
LOOP DROP ;

#64 1- =: szmsk PRIVATE
szmsk 1+ VALUE Z#SP
szmsk 1+ Z#BLOCK 'data PRIVATE

: Z#[] ( addr1 ix -- addr2 ) /MPC * + ;
: []Z# ( ix addr1 -- addr2 ) SWAP Z#[] ;

: Z#PUSH ( -- n ) Z#SP 1- szmsk AND DUP TO Z#SP 'data []Z# ;
: _Z#PICK ( u -- n ) Z#SP + szmsk AND 'data []Z# ; PRIVATE
: Z#DROP ( -- ) Z#SP 1+ szmsk AND TO Z#SP ;
: Z#2DROP ( -- ) Z#SP 2+ szmsk AND TO Z#SP ;
: Z#TOS ( -- n ) 0 _Z#PICK ;
: Z#NOS ( -- n ) 1 _Z#PICK ;
: Z#SEC ( -- n ) 2 _Z#PICK ;
: Z#NSEC ( -- n ) 3 _Z#PICK ;
: Z#TRI ( -- n ) 4 _Z#PICK ;
: Z#NTRI ( -- n ) 5 _Z#PICK ;

-- Input in Z#BASE
: Z#IN ( c-addr u -- )
inumpad SWAP DUP >S CMOVE
S> inumpad + C0!
Z#PUSH inumpad Z#BASE Z#RND mpc_set_str
0< ABORT" Z#IN :: invalid syntax or bad range" ;

: Z#@ ( src -- ) Z#PUSH SWAP Z#RND mpc_set DROP ;
: Z#! ( dst -- ) Z#TOS Z#RND mpc_set DROP Z#DROP ;
: Z#+! ( src -- ) Z#TOS DUP ROT Z#RND mpc_add DROP ;
: Z#-! ( src -- ) Z#TOS DUP ROT Z#RND mpc_sub DROP ;
: Z#*! ( src -- ) Z#TOS DUP ROT Z#RND mpc_mul DROP ;
: Z#/! ( src -- ) Z#TOS DUP ROT Z#RND mpc_div DROP ;

: Z#[@] ( addr1 ix -- ) Z#[] Z#@ ;
: Z#[!] ( addr1 ix -- ) Z#[] Z#! ;
: [@]Z# ( ix addr1 -- ) []Z# Z#@ ;
: [!]Z# ( ix addr1 -- ) []Z# Z#! ;

: Z#DUP ( -- ) Z#TOS Z#@ ;
: Z#OVER ( -- ) Z#NOS Z#@ ;
: Z#SWAP ( -- ) Z#TOS Z#NOS mpc_swap ;
: Z#ROT ( -- ) Z#TOS Z#NOS mpc_swap Z#TOS Z#SEC mpc_swap ;
: -Z#ROT ( -- ) Z#TOS Z#SEC mpc_swap Z#TOS Z#NOS mpc_swap ;
: Z#NIP ( -- ) Z#SWAP Z#DROP ;
: Z#TUCK ( -- ) Z#DUP -Z#ROT ;
: Z#PICK ( n -- ) _Z#PICK Z#@ ;

: Z#2DUP ( -- ) Z#OVER Z#OVER ;
: Z#2OVER ( -- ) 3 Z#PICK 3 Z#PICK ;
: Z#2SWAP ( -- ) Z#NSEC Z#NOS mpc_swap Z#SEC Z#TOS mpc_swap ;
: Z#2NIP ( -- ) Z#2SWAP Z#2DROP ;
: Z#2ROT ( -- ) Z#2SWAP Z#NTRI Z#NOS mpc_swap Z#TRI Z#TOS mpc_swap ;

: Z#$, ( c-addr u -- ) HERE DUP >S /MPC ALLOT Z#PRECISION mpc_init2 Z#IN S> Z#! ;
: Z#VARIABLE ( -- ) 1 Z#BLOCK ;

: Z#CONSTANT ( z: a -- )
CREATE HERE DUP >S /MPC ALLOT
Z#PRECISION mpc_init2 S> Z#!
DOES> Z#@ ;

S" (0 0)" Z#IN Z#CONSTANT 0+0i
S" (0 1)" Z#IN Z#CONSTANT 0+1i
S" (0 -1)" Z#IN Z#CONSTANT 0-1i
S" (1 0)" Z#IN Z#CONSTANT 1+0i
S" (1 1)" Z#IN Z#CONSTANT 1+1i
S" (1 -1)" Z#IN Z#CONSTANT 1-1i

: Z#0! ( dst -- ) 0+0i Z#! ;
: Z#CLEAR ( dst -- ) Z#0! ;

: Z#VALUE ( z: z "str" -- )
CREATE IMMEDIATE
HERE DUP >S /MPC ALLOT Z#PRECISION mpc_init2 S> Z#!
DOES> %VAR @ 0 %VAR !
CASE
( -to) -2 OF ALITERAL EVAL" Z#-!" ENDOF
( +to) -1 OF ALITERAL EVAL" Z#+!" ENDOF
( from) 0 OF ALITERAL EVAL" Z#@" ENDOF
( to) 1 OF ALITERAL EVAL" Z#!" ENDOF
( 0to) 2 OF ALITERAL EVAL" Z#0!" ENDOF
( 'of) 3 OF ALITERAL ENDOF
( sizeof) 4 OF DROP /MPC ILITERAL ENDOF
( /of) 5 OF DROP 1 ILITERAL ENDOF
DUP #-21 ?ERROR" Z#VALUE: undefined message"
ENDCASE ;

: Z#+ ( -- ) Z#NOS DUP Z#TOS Z#RND mpc_add DROP Z#DROP ;
: Z#- ( -- ) Z#NOS DUP Z#TOS Z#RND mpc_sub DROP Z#DROP ;
: Z#* ( -- ) Z#NOS DUP Z#TOS Z#RND mpc_mul DROP Z#DROP ;
: Z#/ ( -- ) Z#NOS DUP Z#TOS Z#RND mpc_div DROP Z#DROP ;

: Z#^N ( n -- ) Z#TOS DUP ROT Z#RND mpc_pow_si DROP ; ( tos <= tos^n )
: Z#POW ( -- ) Z#NOS DUP Z#TOS Z#RND mpc_pow DROP Z#DROP ; ( tos <= nos^tos )
: Z#** ( -- ) Z#POW ; ( tos <= nos^tos )

: Z#2* ( -- ) Z#TOS DUP 2 Z#RND mpc_mul_si DROP ; ( tos <= tos*2 )
: Z#2/ ( -- ) Z#TOS DUP 2 Z#RND mpc_div_ui DROP ; ( tos <= tos/2 )
: Z#N* ( n -- ) Z#TOS DUP ROT Z#RND mpc_mul_si DROP ; ( tos <= tos*n )
: Z#U/ ( u -- ) Z#TOS DUP ROT Z#RND mpc_div_ui DROP ; ( tos <= tos/u )
: Z#N+ ( n -- ) Z#TOS DUP ROT Z#RND mpc_add_si DROP ; ( tos_r <= tos+n )

: Z#1/Z ( -- ) Z#TOS DUP 1 SWAP Z#RND mpc_ui_div DROP ; ( tos <= 1/tos )
: Z#SQR ( -- ) Z#TOS DUP Z#RND mpc_sqr DROP ; ( tos <= sqr[tos] )
: Z#SQRT ( -- ) Z#TOS DUP Z#RND mpc_sqrt DROP ; ( tos <= sqrt[tos] )
: Z#NEGATE ( -- ) Z#TOS DUP Z#RND mpc_neg DROP ; ( tos <= neg[tos] )
: Z#CONJUGATE ( -- ) Z#TOS DUP Z#RND mpc_conj DROP ; ( tos <= conjg[tos] )

: Z#N/ ( n -- ) DUP 0< IF ABS Z#U/ Z#NEGATE EXIT ENDIF Z#U/ ; ( tos <= tos/u )

: Z#ABS ( -- ) F#PUSH Z#TOS Z#RND mpc_abs DROP Z#DROP ; ( ftos <= abs[tos] )
: Z#ABS^2 ( -- ) F#PUSH Z#TOS Z#RND mpc_norm DROP Z#DROP ; ( ftos <= abs^2[tos] )

: Z#SCALE ( -- ) Z#TOS DUP F#TOS Z#RND mpc_mul_fr DROP F#DROP ;
: Z#*F ( -- ) Z#SCALE ;
: Z#+F ( -- ) Z#TOS DUP F#TOS Z#RND mpc_add_fr DROP F#DROP ;
: Z#-F ( -- ) Z#TOS DUP F#TOS Z#RND mpc_sub_fr DROP F#DROP ;
: Z#/F ( -- ) Z#TOS DUP F#TOS Z#RND mpc_div_fr DROP F#DROP ;
: F#/Z ( -- ) Z#TOS F#TOS OVER Z#RND mpc_fr_div DROP F#DROP ;

: Z#EXP(jX) ( -- ) F#SINCOS F#NOS F#TOS Z#PUSH Z#RND mpc_set_fr_fr DROP F#2DROP ;
: Z#EXP(-jX) ( -- ) F#NEGATE Z#EXP(jX) ;

: Z#REAL ( -- ) Z#TOS mpc_realref F#@ Z#DROP ;
: Z#IMAG ( -- ) Z#TOS mpc_imagref F#@ Z#DROP ;

: Z#= ( -- re im ) Z#NOS Z#TOS mpc_cmp DUP 3 AND SWAP 2 RSHIFT Z#2DROP ; \ %iirr; 0 0 is equal

: Z#LN ( -- ) Z#TOS DUP Z#RND mpc_log DROP ; ( tos <= ln[tos] )
: Z#EXP ( -- ) Z#TOS DUP Z#RND mpc_exp DROP ; ( tos <= e^x[tos] )
: Z#SIN ( -- ) Z#TOS DUP Z#RND mpc_sin DROP ; ( tos <= sin[tos] )
: Z#COS ( -- ) Z#TOS DUP Z#RND mpc_cos DROP ; ( tos <= cos[tos] )
: Z#TAN ( -- ) Z#TOS DUP Z#RND mpc_tan DROP ; ( tos <= tan[tos] )
: Z#SINH ( -- ) Z#TOS DUP Z#RND mpc_sinh DROP ; ( tos <= sinh[tos] )
: Z#COSH ( -- ) Z#TOS DUP Z#RND mpc_cosh DROP ; ( tos <= cosh[tos] )
: Z#TANH ( -- ) Z#TOS DUP Z#RND mpc_tanh DROP ; ( tos <= tanh[tos] )
: Z#ASIN ( -- ) Z#TOS DUP Z#RND mpc_asin DROP ; ( tos <= asin[tos] )
: Z#ACOS ( -- ) Z#TOS DUP Z#RND mpc_acos DROP ; ( tos <= acos[tos] )
: Z#ATAN ( -- ) Z#TOS DUP Z#RND mpc_atan DROP ; ( tos <= atan[tos] )
: Z#ASINH ( -- ) Z#TOS DUP Z#RND mpc_asinh DROP ; ( tos <= asinh[tos] )
: Z#ACOSH ( -- ) Z#TOS DUP Z#RND mpc_acosh DROP ; ( tos <= acosh[tos] )
: Z#ATANH ( -- ) Z#TOS DUP Z#RND mpc_atanh DROP ; ( tos <= atanh[tos] )

: Z#SINCOS ( -- ) Z#TOS Z#DUP OVER Z#RND DUP mpc_sin_cos DROP ; ( nos <= sin[tos] tos <= cos[tos] )
: R,I->Z# ( -- ) Z#PUSH DUP mpc_realref SWAP mpc_imagref F#! F#! ;
: Z#->R,I ( -- ) Z#TOS DUP mpc_realref F#@ mpc_imagref F#@ Z#DROP ;

: X+ ( -- ) Z#TOS mpc_realref DUP F#TOS Z#RND mpc_add_fr DROP F#DROP ;
: Y+ ( -- ) Z#TOS mpc_imagref DUP F#TOS Z#RND mpc_add_fr DROP F#DROP ;
: I* ( -- ) Z#->R,I F#NEGATE F#SWAP R,I->Z# ; \ 0+1i Z#*
: CONJ-I* ( -- ) Z#->R,I F#SWAP R,I->Z# ;

: Z#OUT ( -- ) Z#->R,I F#SWAP F#OUT SPACE F#OUT ." i" ;
: Z#OUT.R ( n -- ) #DECIMALS >S TO #DECIMALS Z#OUT S> TO #DECIMALS ;

: .S-PREFIX ( -- ) CR '[' EMIT Z#SP BASE @ >S DECIMAL 2 .R ." ] " S> BASE ! ; PRIVATE

: Z#.S ( -- )
Z#SP 0= IF CR ." empty Z-stack" EXIT ENDIF
Z#SP LOCAL oldsp
BEGIN Z#SP WHILE .S-PREFIX Z#OUT REPEAT
oldsp TO Z#SP ;

: .ZARR ( addr cols -- ) ." [" ( n) 0 ?DO DUP I Z#[@] CR Z#OUT LOOP DROP CR ." ]" ;
: .ZMAT ( addr rows cols -- ) CR SWAP 0 ?DO 2DUP .ZARR TUCK Z#[] SWAP LOOP 2DROP ;

: Z^^ CR Z#.S CR EKEY ^Z = ABORT" huh" ;

:ABOUT CR ." MPC -- arbitrary precision complex FP package based on GMP, MPFR and MPC" ;

NESTING @ 1 = [IF] .ABOUT -mpc CR [THEN]
DEPRIVE

(* End of Source *)