feat(padic): add p-adic numbers

check/README.md

@@ -210,6 +210,8 @@ The following methods are available in Ruby test programs.
   `true` if FORM is the multithreaded version (TFORM), otherwise `false`.
 - `mpi? → bool`  
   `true` if FORM is the MPI version (ParFORM), otherwise `false`.
+- `flint? → bool`  
+  `true` if FORM was built with FLINT support, otherwise `false`.
 - `valgrind? → bool`  
   `true` if FORM is running under Valgrind, otherwise `false`.
 - `wordsize → integer`  

check/check.rb

@@ -310,6 +310,10 @@ def mpi?
     FormTest.cfg.mpi?
   end
 
+  def flint?
+    FormTest.cfg.flint?
+  end
+
   def valgrind?
     if FormTest.cfg.fake_valgrind.nil?
       !FormTest.cfg.valgrind.nil?
@@ -1249,6 +1253,7 @@ def initialize(form, mpirun, mpirun_opts, valgrind, valgrind_opts, fake_valgrind
     @is_serial   = nil
     @is_threaded = nil
     @is_mpi      = nil
+    @has_flint   = nil
     @wordsize    = wordsize
     @form_cmd    = nil
   end
@@ -1268,6 +1273,10 @@ def mpi?
     @is_mpi
   end
 
+  def flint?
+    @has_flint
+  end
+
   def check_bin(name, bin)
     # Check if the executable is available.
     if File.executable?(bin)
@@ -1326,6 +1335,17 @@ def check
         system("#{form_bin} #{frmname}")
         fatal("failed to get the version of '#{@form}'")
       end
+      # Check if Flint is available.
+      # Use form -vv which contains either "+flint=VERSION" or "-flint".
+      feature_out, _status = Open3.capture2e(@form_bin, "-vv")
+      if feature_out =~ /(?:^|\s)\+flint(?:=|\s|$)/
+        @has_flint = true
+      elsif feature_out =~ /(?:^|\s)-flint(?:\s|$)/
+        @has_flint = false
+      else
+        warn("failed to determine FLINT support of '#{@form}'")
+        @has_flint = false
+      end
       # Check the wordsize.
       # Method 1: from the output header.
       # Method 2: 2^64 = 0 (mod 2^64) => 64-bit machine => sizeof(WORD) = 4, etc.

check/features.frm

@@ -2382,6 +2382,157 @@ Evaluate euler_;
 #pend_if wordsize == 2
 assert runtime_error?("Divergent Euler sum in CalculateEuler")
 *--#] mzv_error_5 :
+*--#[ padic_1 :
+#StartPadic 7,10
+Local F = 101/7;
+Local G = 7/13;
+ToPadic;
+.sort
+
+Local MUL = (10+F)*G;
+Local DIV = (10+F)/G;
+Print;
+.end
+#require flint?
+#pend_if wordsize == 2
+assert succeeded?
+assert result("F") =~ expr("(3*7^-1 + 2*7^1)")
+assert result("G") =~ expr("(6*7^1 + 4*7^2 + 2*7^3 + 5*7^4 + 3*7^5 + 1*7^7 + 2*7^8 + 4*7^9)")
+assert result("MUL") =~ expr("(4 + 4*7^1 + 5*7^2 + 3*7^3 + 1*7^5 + 2*7^6 + 4*7^7 + 1*7^8)")
+assert result("DIV") =~ expr("(4*7^-2 + 2*7^-1 + 3 + 6*7^1 + 3*7^7 + 1*7^9)")
+*--#] padic_1 :
+*--#[ padic_2 :
+#-
+Off Statistics;
+#StartPadic 7,10
+Symbol a,x1,...,x4;
+Function A1,A2;
+CFunction f1,f2;
+Vector p1,p2,q;
+Local F = 31*x1*x2*f2*f1*q*p1.p2*A1*A2*cos_(3*pi_)*mzv_(2,1)*mzv_*cos_*mzv_(2)*cos_(pi_);
+ToPadic;
+.sort
+
+* Notice that the padic_ function is the last element in the term. 
+* See normalize.c for the order of the terms:
+* 1. Non-commuting functions
+* 2. Commuting functions
+* 3. LEVICIVITA tensors.
+* 4. DELTA.
+* 5. loose INDEX.
+* 6. VECTOR.
+* 7. DOTPRODUCT.
+* 8. SYMBOL.
+* 9. optional float_ / padic_ coefficient functions, if enabled.
+* 10. rational coefficient at the very end.
+* Notice that functions are first ordered according to their function number, 
+* see ftypes.h for the predefined functions, and then according to their arguments. 
+Format PadicPrint off;
+Print;
+.sort
+
+#EndPadic
+* padic_ should now not be placed at the very end of the term, 
+* but together with the other non-commuting functions.
+* At the moment, padic_ is the last predefined function, but it's number
+* is of course lower than that of the first user defined function. 
+Local G = F;
+Print G;
+.end
+#require flint?
+#pend_if wordsize == 2
+assert succeeded?
+assert result("F") =~ expr("A1*A2*cos_*cos_(pi_)*cos_(3*pi_)*mzv_*mzv_(2)*mzv_(2,1)*f1*f2*q*p1.p2*x1
+      *x2*padic_(0,10,31)")
+assert result("G") =~ expr("A1*A2*cos_*cos_(pi_)*cos_(3*pi_)*mzv_*mzv_(2)*mzv_(2,1)*padic_(0,10,31)*
+      f1*f2*q*p1.p2*x1*x2")
+*--#] padic_2 :
+*--#[ AddWithPadic :
+#-
+Off Statistics;
+#StartPadic 7,10
+Symbol x1,...,x4;
+
+Local F = (x1+padic_(0,10,1)*x2+x3+padic_(0,10,1)*x4)^5;
+id x1 = 1-x2-x3-x4;
+Print;
+.end
+#require flint?
+#pend_if wordsize == 2
+assert succeeded?
+assert result("F") =~ expr("1")
+*--#] AddWithPadic :
+*--#[ MergeWithPadic :
+#: termsinsmall 16
+On fewerstats 1;
+#StartPadic 7,10
+Format PadicPrint off;
+Symbol x1,...,x16;
+
+Local F = (x1+padic_(0,10,1)*x2+x3+padic_(0,10,1)*x4)^5;
+id x1 = 1-x2-x3-x4;
+.sort
+
+Local CORNER = x1+...+x15+padic_(0,10,1)*x16+2^80*x16;
+Print;
+.end
+#require flint?
+#pend_if wordsize == 2
+assert succeeded?
+assert result("F") =~ expr("1")
+assert result("CORNER") =~ expr("x16*padic_(0,10,188220513) + x15 + x14 + x13 + x12 + x11 + x10 + x9 + x8
+       + x7 + x6 + x5 + x4 + x3 + x2 + x1")
+*--#] MergeWithPadic :
+*--#[ padic_reconstruct_1 :
+Symbol x,j;
+#StartPadic 7,5
+Local F = 101/13*sum_(j,-5,5,x^j*7^j);
+ToPadic;
+.sort
+
+PadicToRat;
+Print;
+.end
+#require flint?
+#pend_if wordsize == 2
+assert succeeded?
+assert result("F") =~ expr("101/13 + 101/218491*x^-5 + 101/31213*x^-4 + 101/4459*x^-3 + 101/637*x^-2
+       + 101/91*x^-1 - 42/23*x - 833/8*x^2 + 1372*x^3 + 2401/2*x^4")
+*--#] padic_reconstruct_1 :
+*--#[ padic_reconstruct_2 :
+#-
+#$prime = 175519;
+#define EP "`$prime'"
+
+Off Statistics;
+Symbol eps,x1,x2,v,N;
+CFunction padic,rat,coeff;
+Table ep();
+Fill ep = `EP';
+
+#StartPadic `$prime',5
+Local F = rat(70944*ep^5 - 20088*ep^4 - 5580*ep^3 + 1489*ep^2 + 126
+      *ep - 31,32000*ep^5 - 3280*ep^3 + 80*ep);
+id rat(x1?,x2?) = x1/x2;
+ToPadic;
+.sort
+
+#do i=0,5
+    FromPadic, padic;
+    id padic(?x1) = padic(?x1)+padic_(?x1);
+    Transform decode(13,3):base=`$prime';
+    id padic(v?,N?,x1?,?x2) = coeff(v,x1)-makerational_(x1,$prime)*`EP'^v;
+    .sort
+#enddo
+id coeff(v?,x1?) = makerational_(x1,$prime)*eps^v;
+.sort
+Print;
+.end
+#require flint?
+#pend_if wordsize == 2
+assert succeeded?
+assert result("F") =~ expr("63/40 - 31/80*eps^-1 + 109/40*eps - 207/40*eps^2 + 125/8*eps^3 + 357/8*eps^4")
+*--#] padic_reconstruct_2 : 
 *--#[ humanstats :
 #-
 On humanstats;

configure.ac

@@ -387,7 +387,7 @@ AS_IF([test -f "$srcdir/extern/zstd/zlibWrapper/zstd_zlibwrapper.h"], [],
 # Check for flint
 AC_ARG_WITH([flint],
 	[AS_HELP_STRING([--with-flint@<:@=DIR@:>@],
-		[use flint for polynomial arithmetic (installed in prefix DIR) @<:@default=check@:>@])],
+		[use flint for polynomial arithmetic and p-adics (installed in prefix DIR) @<:@default=check@:>@])],
 	[AS_IF([test "x$withval" != xyes && test "x$withval" != xno && test "x$withval" != xcheck],
 		[with_flint=yes
 		CPPFLAGS="$CPPFLAGS -I$withval/include"
@@ -417,7 +417,7 @@ AS_IF([test "x$with_flint" != xno],
 				[MATH_LIBS="$ac_cv_search_pow"])])])
 	AS_IF([$flag], [AC_CHECK_LIB([flint], [fmpz_mpoly_init], [LIBS="-lflint $MATH_LIBS $LIBS"], [flag=false], [$MATH_LIBS])])
 	AS_IF([$flag],
-		[AC_DEFINE(WITHFLINT, [], [Define to use flint for polynomial arithmetic.])
+		[AC_DEFINE(WITHFLINT, [], [Define to use flint for polynomial arithmetic and p-adics.])
 		with_flint=yes],
 		[AS_IF([test "x$with_flint" = xyes],
 			[AC_MSG_FAILURE([test for flint failed. Give --without-flint if you want to compile without flint])])

sources/Makefile.am

@@ -89,7 +89,8 @@ if WITHFLINT
 SRCBASE += \
 	flintinterface.cc \
 	flintinterface.h \
-	flintwrap.cc
+	flintwrap.cc \
+	padic.c
 endif
 
 if WITHFLOAT

sources/argument.c

@@ -171,6 +171,9 @@ int execarg(PHEAD WORD *term, WORD level)
 					if ( m[1] == 2 ) {
 #ifdef WITHFLOAT
 						if ( *t != FLOATFUN || TestFloat(t) == 0 )
+#endif
+#ifdef WITHFLINT
+						if ( *t != PADICFUN || TestPadic(t) == 0 )
 #endif
 						{
 							m += 2;
@@ -181,6 +184,11 @@ int execarg(PHEAD WORD *term, WORD level)
 						else {
 							m += 2;
 						}
+#endif
+#ifdef WITHFLINT
+						else {
+							m += 2;
+						}
 #endif
 					}
 					else {

sources/compcomm.c

@@ -58,6 +58,9 @@ static KEYWORD formatoptions[] = {
 	,{"mathematica",	(TFUN)0,	MATHEMATICAMODE,	0}
 	,{"normal",			(TFUN)0,	NORMALFORMAT,		1}
 	,{"nospaces",		(TFUN)0,	NOSPACEFORMAT,		3}
+#ifdef WITHFLINT
+	,{"padicprint",		(TFUN)0,	0,					6}
+#endif
 	,{"pfortran",		(TFUN)0,	PFORTRANMODE,		0}
 	,{"quadfortran",	(TFUN)0,	QUADRUPLEFORTRANMODE,	0}
 	,{"quadruplefortran",	(TFUN)0,	QUADRUPLEFORTRANMODE,	0}
@@ -410,6 +413,39 @@ WrongOption:		MesPrint("&Illegal option in Format FloatPrecision: %s",s);
 					error = 1;
 				}
 			}
+#endif
+#ifdef WITHFLINT
+			else if ( key->flags == 6 ) {
+/*
+				Syntax: Format PadicPrint [on];
+				        Format PadicPrint off;
+*/
+				while ( FG.cTable[*s] == 0 ) s++;
+				while ( *s == ' ' || *s == '\t' || *s == ',' ) s++;
+				if ( *s == 0 ) {
+					AO.PadicPrint = 1;
+				}
+				else if ( tolower(*s) == 'o' && tolower(s[1]) == 'n' ) {
+					ss = s;
+					s += 2;
+					while ( *s == ' ' || *s == '\t' || *s == ',' ) s++;
+					if ( *s ) { s = ss; goto WrongPadicOption; }
+					AO.PadicPrint = 1;
+				}
+				else if ( tolower(*s) == 'o' && tolower(s[1]) == 'f'
+				&& tolower(s[2]) == 'f' ) {
+					ss = s;
+					s += 3;
+					while ( *s == ' ' || *s == '\t' || *s == ',' ) s++;
+					if ( *s ) { s = ss; goto WrongPadicOption; }
+					AO.PadicPrint = 0;
+				}
+				else {
+WrongPadicOption:
+					MesPrint("&Illegal option in Format PadicPrint: %s",s);
+					error = 1;
+				}
+			}
 #endif
 		}
 		else if ( ( *s == 'c' || *s == 'C' ) && ( FG.cTable[s[1]] == 1 ) ) {
@@ -2067,6 +2103,19 @@ doset:					if ( Sets[number].type != CFUNCTION ) goto nofun;
 								error = 1;
 							}
 						}
+#endif
+#ifdef WITHFLINT
+						{
+							WORD *r1, *r2;
+							r1 = SetElements + Sets[number].first;
+							r2 = SetElements + Sets[number].last;
+							while ( r1 < r2 ) {
+								if ( *r1++ == PADICFUN ) {
+									MesPrint("&Illegal use of argument environment and padic_.");
+									error = 1;
+								}
+							}
+						}
 #endif
 						*w++ = CSET; *w++ = number;
 					}
@@ -2076,6 +2125,12 @@ doset:					if ( Sets[number].type != CFUNCTION ) goto nofun;
 							MesPrint("&Illegal use of argument environment and float_.");
 							error = 1;
 						}
+#endif
+#ifdef WITHFLINT
+						if ( (number + FUNCTION) == PADICFUN ) {
+							MesPrint("&Illegal use of argument environment and padic_.");
+							error = 1;
+						}
 #endif
 						*w++ = CFUNCTION; *w++ = number + FUNCTION;
 					}
@@ -3599,6 +3654,12 @@ int CoModulus(UBYTE *inp)
 		MesPrint("&Simultaneous use of floating point and modulus arithmetic makes no sense.");
 		Retval = 1;
 	}
+#endif
+#ifdef WITHFLINT
+	if ( AC.activePadic ) {
+		MesPrint("&Simultaneous use of p-adic and modulus arithmetic makes no sense.");
+		Retval = 1;
+	}
 #endif
 	AC.modmode = 0;
 	if ( *inp == '-' ) {
@@ -5388,6 +5449,12 @@ int CoPolyFun(UBYTE *s)
 		MesPrint("&Simultaneous use of PolyFun and float_ is not allowed.");
 		error = 1;
 	}
+#endif
+#ifdef WITHFLINT
+	if ( AC.activePadic ) {
+		MesPrint("&Simultaneous use of PolyFun and padic_ is not allowed.");
+		error = 1;
+	}
 #endif
 	return(error);
 }
@@ -5420,6 +5487,12 @@ int CoPolyRatFun(UBYTE *s)
 		MesPrint("&Simultaneous use of PolyFun and float_ is not allowed.");
 		error = 1;
 	}
+#endif
+#ifdef WITHFLINT
+	if ( AC.activePadic ) {
+		MesPrint("&Simultaneous use of PolyFun and padic_ is not allowed.");
+		error = 1;
+	}
 #endif
 	if ( ( ( type = GetName(AC.varnames,s,&numfun,WITHAUTO) ) != CFUNCTION )
 	|| ( functions[numfun].spec != 0 ) || ( functions[numfun].commute != 0 ) ) {