Add MPFR support to ForeignFunctions package

M2/Macaulay2/d/ffi.d

@@ -207,13 +207,18 @@ setupconst("ffiVoidType", toExpr(Ccode(voidPointer, "&ffi_type_void")));
 -- integer types --
 -------------------
 
+-- returns pointer to ffi_type object for given integer type
+-- input: (n:ZZ, signed:Boolean)
+-- n = number of bits (or 0 for mpz_t)
 ffiIntegerType(e:Expr):Expr := (
     when e
     is a:Sequence do (
 	if length(a) == 2 then (
 	    when a.0
 	    is n:ZZcell do (
-		when a.1
+		if isZero(n.v)
+		then toExpr(Ccode(voidPointer, "&ffi_type_pointer"))
+		else when a.1
 		is signed:Boolean do (
 		    bits := toInt(n);
 		    if signed.v then (
@@ -243,23 +248,24 @@ ffiIntegerType(e:Expr):Expr := (
     else WrongNumArgs(2));
 setupfun("ffiIntegerType", ffiIntegerType);
 
+-- returns pointer to integer object with given value
+-- inputs (x:ZZ, y:ZZ, signed:Boolean)
+-- x = value
+-- y = number of bits (or 0 for mpz_t)
 ffiIntegerAddress(e:Expr):Expr := (
     when e
-    is x:ZZcell do (
-	y := newZZmutable();
-	set(y, x.v);
-	ptr := getMem(pointerSize);
-	Ccode(void, "*(void **)", ptr, " = ", y);
-	Ccode(void, "GC_REGISTER_FINALIZER(", ptr,
-	    ", (GC_finalization_proc)mpz_clear, ", y, ", 0, 0)");
-	toExpr(ptr))
     is a:Sequence do (
 	if length(a) == 3 then (
 	    when a.0
 	    is x:ZZcell do (
 		when a.1
 		is y:ZZcell do (
-		    when a.2
+		    if isZero(y.v) then (
+			z := copy(x.v);
+			ptr := getMem(pointerSize);
+			Ccode(void, "*(mpz_srcptr *)", ptr, " = ", z);
+			toExpr(ptr))
+		    else when a.2
 		    is signed:Boolean do (
 			bits := toInt(y);
 			ptr := getMemAtomic(bits / 8);
@@ -302,16 +308,21 @@ ffiIntegerAddress(e:Expr):Expr := (
     else WrongNumArgs(3));
 setupfun("ffiIntegerAddress", ffiIntegerAddress);
 
+-- returns value (as ZZ) of integer object given its address
+-- inputs: (x:Pointer, y:ZZ, signed:ZZ)
+-- x = pointer to integer object
+-- y = number of bits (or 0 for mpz_t)
 ffiIntegerValue(e:Expr):Expr := (
     when e
-    is x:pointerCell do toExpr(moveToZZ(Ccode(ZZmutable, "*(mpz_ptr*)", x.v)))
     is a:Sequence do (
 	if length(a) == 3 then (
 	    when a.0
 	    is x:pointerCell do (
 		when a.1
 		is y:ZZcell do (
-		    when a.2
+		    if isZero(y.v)
+		    then toExpr(moveToZZ(Ccode(ZZmutable, "*(mpz_ptr*)", x.v)))
+		    else when a.2
 		    is signed:Boolean do (
 			bits := toInt(y);
 			if signed.v then (
@@ -347,16 +358,23 @@ setupfun("ffiIntegerValue", ffiIntegerValue);
 -- real types --
 ----------------
 
+-- returns pointer to ffi_type object for given real number type
+-- input: n:ZZ (0 = mpfr_t, 32 = float, 64 = double)
 ffiRealType(e:Expr):Expr := (
     when e
     is n:ZZcell do (
 	bits := toInt(n);
-	if bits == 32 then toExpr(Ccode(voidPointer, "&ffi_type_float"))
+	if bits == 0 then toExpr(Ccode(voidPointer, "&ffi_type_pointer"))
+	else if bits == 32 then toExpr(Ccode(voidPointer, "&ffi_type_float"))
 	else if bits == 64 then toExpr(Ccode(voidPointer, "&ffi_type_double"))
 	else buildErrorPacket("expected 32 or 64 bits"))
     else WrongArgZZ());
 setupfun("ffiRealType", ffiRealType);
 
+-- returns pointer to real number object with given value
+-- inputs: (x:RR, y:ZZ)
+-- x = value
+-- y = type (0 = mpfr_t, 32 = float, 64 = double)
 ffiRealAddress(e:Expr):Expr := (
     when e
     is a:Sequence do (
@@ -366,23 +384,33 @@ ffiRealAddress(e:Expr):Expr := (
 		when a.1
 		is y:ZZcell do (
 		    bits := toInt(y);
-		    ptr := getMemAtomic(bits / 8);
-		    if bits == 32
-		    then (
-			z := toFloat(x);
-			Ccode(void, "*(float *)", ptr, " = ", z))
-		    else if bits == 64
-		    then (
-			z := toDouble(x);
-			Ccode(void, "*(double *)", ptr, " = ", z))
-		    else return buildErrorPacket("expected 32 or 64 bits");
-		    toExpr(ptr))
+		    if bits == 0 then (
+			z := copy(x.v);
+			ptr := getMem(pointerSize);
+			Ccode(void, "*(mpfr_srcptr *)", ptr, " = ", z);
+			toExpr(ptr))
+		    else if bits == 32 || bits == 64 then (
+			ptr := getMemAtomic(bits / 8);
+			if bits == 32
+			then (
+			    z := toFloat(x);
+			    Ccode(void, "*(float *)", ptr, " = ", z))
+			else if bits == 64
+			then (
+			    z := toDouble(x);
+			    Ccode(void, "*(double *)", ptr, " = ", z));
+			toExpr(ptr))
+		    else return buildErrorPacket("expected 0, 32, or 64"))
 		else WrongArgZZ(2))
 	    else WrongArgRR(1))
 	else WrongNumArgs(2))
     else WrongNumArgs(2));
 setupfun("ffiRealAddress", ffiRealAddress);
 
+-- returns value (as RR) of real number object given its address
+-- inputs: (x:Pointer, y:ZZ)
+-- x = pointer to real number object
+-- y = type (0 = mpfr_t, 32 = float, 64 = double)
 ffiRealValue(e:Expr):Expr := (
     when e
     is a:Sequence do (
@@ -392,11 +420,13 @@ ffiRealValue(e:Expr):Expr := (
 		when a.1
 		is y:ZZcell do (
 		    bits := toInt(y);
-		    if bits == 32
+		    if bits == 0
+		    then toExpr(moveToRR(Ccode(RRmutable, "*(mpfr_ptr*)", x.v)))
+		    else if bits == 32
 		    then toExpr(Ccode(float, "*(float *)", x.v))
 		    else if bits == 64
 		    then toExpr(Ccode(double, "*(double *)", x.v))
-		    else buildErrorPacket("expected 32 or 64 bits"))
+		    else buildErrorPacket("expected 0, 32, or 64"))
 		else WrongArgZZ(2))
 	    else WrongArgPointer(1))
 	else WrongNumArgs(2))

M2/Macaulay2/d/gmp.d

@@ -203,7 +203,7 @@ clear(x:RRimutable) ::= Ccode( void, "mpfi_clear(",  x, ")" );
 
 clear(z:CCmutable):void := ( clear(z.re); clear(z.im); );
 
-export moveToZZ(z:ZZmutable):ZZ := (
+export copy(z:ZZ):ZZ := (
      y := GCmalloc(ZZmutable);
      Ccode(void, "
 	  int s = z->_mp_size, ss = s>=0 ? s : -s;
@@ -212,13 +212,14 @@ export moveToZZ(z:ZZmutable):ZZ := (
 	  ",y,"->_mp_alloc = ss, ",y,"->_mp_size = s, ",y,"->_mp_d = p;
 	  ");
      Ccode(ZZ,y));
+export moveToZZ(z:ZZmutable):ZZ := copy(Ccode(ZZ, z));
 
 export moveToZZandclear(z:ZZmutable):ZZ := (
      w := moveToZZ(z);
      clear(z);
      w);
 
-export moveToRR(z:RRmutable):RR := (
+export copy(z:RR):RR := (
      y := GCmalloc(RRmutable);
      Ccode(void, "
   	  int limb_size = (",z,"->_mpfr_prec - 1) / GMP_NUMB_BITS + 1;
@@ -231,6 +232,7 @@ export moveToRR(z:RRmutable):RR := (
 	  ");
     Ccode(RR,y)
    );
+export moveToRR(z:RRmutable):RR := copy(Ccode(RR, z));
 
 export moveToRRi(z:RRimutable):RRi := (
     y := GCmalloc(RRimutable);

M2/Macaulay2/packages/ForeignFunctions.m2

@@ -1,7 +1,7 @@
 newPackage("ForeignFunctions",
     Headline => "foreign function interface",
     Version => "0.3",
-    Date => "January 28, 2023",
+    Date => "February 8, 2024",
     Authors => {{
 	    Name => "Doug Torrance",
 	    Email => "dtorrance@piedmont.edu",
@@ -15,9 +15,9 @@ newPackage("ForeignFunctions",
 
 -*
 
-0.3 (2024-01-28, M2 1.23)
+0.3 (2024-02-08, M2 1.23)
 * add subscripted assignment for various pointer types
-* add support for GMP integers
+* add support for GMP integers and MPFR reals
 * add support for describe, expression, toExternalString, and toString
 * use null coalescing operator
 
@@ -77,6 +77,7 @@ export {
     "mpzT",
     "float",
     "double",
+    "mpfrT",
     "voidstar",
     "charstar",
     "voidstarstar",
@@ -258,12 +259,7 @@ if version#"pointer size" == 4 then (
     ) else (
     long = int64;
     ulong = uint64)
-
-mpzT = new ForeignIntegerType;
-mpzT.Name = "mpz_t";
-mpzT.Address = ffiPointerType
-new mpzT from ZZ := (T, n) -> new T from {Address => ffiIntegerAddress n}
-value mpzT := ffiIntegerValue @@ address
+mpzT = foreignIntegerType("mpz_t", 0, true)
 
 ForeignIntegerType Number :=
 ForeignIntegerType Constant := (T, x) -> new T from truncate x
@@ -289,12 +285,13 @@ foreignRealType(String, ZZ) := (name, bits) -> (
 
 float = foreignRealType("float", 32)
 double = foreignRealType("double", 64)
+mpfrT = foreignRealType("mpfr_t", 0)
 
 ForeignRealType Number :=
 ForeignRealType Constant := (T, x) -> new T from realPart numeric x
 ForeignRealType RRi := (T, x) -> T toRR x
 
-isAtomic ForeignRealType := T -> true
+isAtomic ForeignRealType := T -> if T === mpfrT then false else true
 
 --------------------------
 -- foreign pointer type --
@@ -957,6 +954,29 @@ doc ///
       double
 ///
 
+doc ///
+  Key
+    mpfrT
+    (NewFromMethod, mpfrT, RR)
+    (value, mpfrT)
+  Headline
+    MPFR multiple-precision floating-point type
+  Description
+    Text
+      Macaulay2's native @TO RR@ real number type wraps around @TT
+      "mpfr_t"@ from @HREF{"https://mpfr.org/", "MPFR"}@.  This type
+      (which is an instance of @TO ForeignRealType@) allows for
+      conversion between Macaulay2 reals and MPFR reals without loss
+      of precision.
+    Example
+      mpfrT numeric(100, pi)
+      value oo
+      mpfrAdd = foreignFunction("mpfr_add", void, {mpfrT, mpfrT, mpfrT, int})
+      x = mpfrT 0p100
+      mpfrAdd(x, numeric(100, pi), exp 1p100, 0)
+      x
+///
+
 doc ///
   Key
     (symbol SPACE, ForeignRealType, Number)
@@ -1962,6 +1982,7 @@ assert Equation(value mpzT 10^100, 10^100)
 -- real types
 assert Equation(value float 3.14159, 3.14159p24)
 assert Equation(value double 3.14159, 3.14159p53)
+assert Equation(value mpfrT 3.14159p100, 3.14159p100)
 
 -- pointer types
 ptr = address int 3