diff --git a/check/README.md b/check/README.md index a9a36c65..268d72e9 100644 --- a/check/README.md +++ b/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` diff --git a/check/check.rb b/check/check.rb index 09d7b32a..835e5407 100755 --- a/check/check.rb +++ b/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. diff --git a/check/features.frm b/check/features.frm index aafc5384..8879aa75 100644 --- a/check/features.frm +++ b/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; diff --git a/configure.ac b/configure.ac index d40185a5..41e2fd3e 100644 --- a/configure.ac +++ b/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])]) diff --git a/sources/Makefile.am b/sources/Makefile.am index d8c3caef..4bda768e 100644 --- a/sources/Makefile.am +++ b/sources/Makefile.am @@ -89,7 +89,8 @@ if WITHFLINT SRCBASE += \ flintinterface.cc \ flintinterface.h \ - flintwrap.cc + flintwrap.cc \ + padic.c endif if WITHFLOAT diff --git a/sources/argument.c b/sources/argument.c index 9e972d5f..f0f7f524 100644 --- a/sources/argument.c +++ b/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 { diff --git a/sources/compcomm.c b/sources/compcomm.c index ac38ead5..ddde7940 100644 --- a/sources/compcomm.c +++ b/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 ) ) { diff --git a/sources/compiler.c b/sources/compiler.c index e7e0fc05..d7a000b6 100644 --- a/sources/compiler.c +++ b/sources/compiler.c @@ -162,6 +162,9 @@ static KEYWORD com2commands[] = { ,{"factorize", (TFUN)CoFactorize, TOOUTPUT, PARTEST} ,{"fill", (TFUN)CoFill, DECLARATION, PARTEST} ,{"fillexpression", (TFUN)CoFillExpression, DECLARATION, PARTEST} +#ifdef WITHFLINT + ,{"frompadic", (TFUN)CoFromPadic, STATEMENT, PARTEST} +#endif ,{"frompolynomial", (TFUN)CoFromPolynomial, STATEMENT, PARTEST} ,{"funpowers", (TFUN)CoFunPowers, DECLARATION, PARTEST} ,{"hide", (TFUN)CoHide, SPECIFICATION,PARTEST} @@ -197,6 +200,9 @@ static KEYWORD com2commands[] = { ,{"on", (TFUN)CoOn, DECLARATION, PARTEST} ,{"once", (TFUN)CoOnce, STATEMENT, PARTEST} ,{"only", (TFUN)CoOnly, STATEMENT, PARTEST} +#ifdef WITHFLINT + ,{"padictorat", (TFUN)CoPadicToRat, STATEMENT, PARTEST} +#endif ,{"particle", (TFUN)CoParticle, DECLARATION, PARTEST} ,{"polyfun", (TFUN)CoPolyFun, DECLARATION, PARTEST} ,{"polyratfun", (TFUN)CoPolyRatFun, DECLARATION, PARTEST} @@ -236,6 +242,9 @@ static KEYWORD com2commands[] = { ,{"threadbucketsize",(TFUN)CoThreadBucket, DECLARATION, PARTEST} #ifdef WITHFLOAT ,{"tofloat", (TFUN)CoToFloat, STATEMENT, PARTEST} +#endif +#ifdef WITHFLINT + ,{"topadic", (TFUN)CoToPadic, STATEMENT, PARTEST} #endif ,{"topolynomial", (TFUN)CoToPolynomial, STATEMENT, PARTEST} #ifdef WITHFLOAT diff --git a/sources/declare.h b/sources/declare.h index 8056b656..9037ee41 100644 --- a/sources/declare.h +++ b/sources/declare.h @@ -366,7 +366,6 @@ TP=T+1;while(TP 0 ) goto NormZero; + } + else { + padicret = MulPadics(BHEAD padicaccum,padicaccum,t); + if ( padicret < 0 ) goto FromNorm; + if ( padicret > 0 ) goto NormZero; + } + withpadic++; + } + break; #endif case INTFUNCTION : /* @@ -2787,6 +2826,34 @@ TryAgain:; stop = m + 3; r = tt; } +#endif +#ifdef WITHFLINT + else if ( *t == PADICFUN && TestPadic(t) ) { + k = t[1]; + pden[i][1] -= k; + pden[i][FUNHEAD] -= k; + pden[i][FUNHEAD+ARGHEAD] -= k; + if ( withpadic == 0 ) { + padicret = InvPadic(BHEAD padicaccum,t); + if ( padicret < 0 ) goto FromNorm; + if ( padicret > 0 ) goto NormZero; + withpadic = 2; + } + else { + padicret = DivPadics(BHEAD padicaccum, + (withpadic == 1) ? firstpadic : padicaccum,t); + if ( padicret < 0 ) goto FromNorm; + if ( padicret > 0 ) goto NormZero; + withpadic++; + } + tt = to = t; + while ( r < m ) *to++ = *r++; + *to++ = 1; *to++ = 1; *to++ = 3; + if ( ncoef < 0 ) t[-1] = -t[-1]; + m -= k; + stop = m + 3; + r = tt; + } #endif t = r; } @@ -4109,6 +4176,28 @@ regularratfun:; } else AT.FloatPos = 0; #endif + #ifdef WITHFLINT + if ( withpadic ) { + WORD *padicfunction = (withpadic == 1) ? firstpadic : padicaccum; + /* + First check whether the coefficient is already 1/1. + */ + if ( ncoef == 3 && n_coef[0] == 1 && n_coef[1] == 1 ) { + AT.PadicPos = m-termout; + i = padicfunction[1]; + NCOPY(m,padicfunction,i) + } + else { + padicret = MulRatToPadic(BHEAD m,padicfunction,(UWORD *)n_coef,ncoef); + if ( padicret < 0 ) goto FromNorm; + if ( padicret > 0 ) goto NormZero; + AT.PadicPos = m-termout; + m += m[1]; + } + n_coef[0] = 1; n_coef[1] = 1; ncoef = 3; + } + else AT.PadicPos = 0; +#endif /* #] float_ : @@ -4163,6 +4252,9 @@ regularratfun:; AT.NormDepth--; TermFree(n_llnum,"n_llnum"); TermFree(n_coef,"NormCoef"); +#ifdef WITHFLINT + TermFree(padicaccum,"NormPadic"); +#endif return(1); } else { @@ -4191,6 +4283,9 @@ regularratfun:; AT.NormDepth--; TermFree(n_llnum,"n_llnum"); TermFree(n_coef,"NormCoef"); +#ifdef WITHFLINT + TermFree(padicaccum,"NormPadic"); +#endif return(regval); NormInf: @@ -4217,12 +4312,18 @@ regularratfun:; AT.NormDepth--; TermFree(n_llnum,"n_llnum"); TermFree(n_coef,"NormCoef"); +#ifdef WITHFLINT + TermFree(padicaccum,"NormPadic"); +#endif return(regval); NormMin: AT.NormDepth--; TermFree(n_llnum,"n_llnum"); TermFree(n_coef,"NormCoef"); +#ifdef WITHFLINT + TermFree(padicaccum,"NormPadic"); +#endif return(-1); FromNorm: @@ -4232,6 +4333,9 @@ regularratfun:; AT.NormDepth--; TermFree(n_llnum,"n_llnum"); TermFree(n_coef,"NormCoef"); +#ifdef WITHFLINT + TermFree(padicaccum,"NormPadic"); +#endif return(-1); /* diff --git a/sources/padic.c b/sources/padic.c new file mode 100644 index 00000000..717aaafe --- /dev/null +++ b/sources/padic.c @@ -0,0 +1,1389 @@ +/** @file padic.c + * + * This file contains the p-adic runtime integration. + * The implementation follows the same integration points as float.c: + * - a dedicated internal function `padic_` used as coefficient carrier, + * - statement support (`ToPadic`) in compiler/executor, + * - normalization/sorting helper routines for coefficient arithmetic, + * - print support with `Format padicprint`. + * + * The numerical backend is FLINT's padic module. + */ +/* #[ License : */ +/* + * Copyright (C) 1984-2026 J.A.M. Vermaseren + * When using this file you are requested to refer to the publication + * J.A.M.Vermaseren "New features of FORM" math-ph/0010025 + * This is considered a matter of courtesy as the development was paid + * for by FOM the Dutch physics granting agency and we would like to + * be able to track its scientific use to convince FOM of its value + * for the community. + * + * This file is part of FORM. + * + * FORM is free software: you can redistribute it and/or modify it under the + * terms of the GNU General Public License as published by the Free Software + * Foundation, either version 3 of the License, or (at your option) any later + * version. + * + * FORM 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 General Public License for more + * details. + * + * You should have received a copy of the GNU General Public License along + * with FORM. If not, see . + */ +/* #] License : */ +/* + #[ Includes: padic.c +*/ +#include "form3.h" +#include +#include + +/* + * FLINT's padic mpq conversion API requires gmp.h to be included before + * flint headers. + */ +#include +#if defined(WINDOWS) +// flint.h defines WORD(xx), which conflicts with the one defined in form3.h. +// see also flintinterface.h +#undef WORD +#endif + +#include +#include +#include +#include + +#if defined(WINDOWS) +// Redefine WORD here to match form3.h. +#undef WORD +#define WORD FORM_WORD +#endif + +/* + FORM keeps a single active p-adic context per run: + - prime p + - precision N + The context is configured by #StartPadic / #EndPadic. + + Note: the internal `padic_` coefficient carrier does not store p itself. + The active `PadicContext` is assumed when unpacking and + operating on p-adic coefficients. +*/ +#define PadicActive AC.activePadic +#define PadicPrecision AC.activePadicPrecision +#define ActivePadicContext AC.activePadicContext +#define PadicContext ((padic_ctx_struct *)(ActivePadicContext)) + +/* + Thread-local aux storage. The pointer is stored in AT.padic_aux_. +*/ +typedef struct PADIC_AUX_ { + /* + All p-adic operations are done through this per-thread scratch space. + This avoids repeated init/clear churn while sorting/normalizing and keeps + GMP/FLINT temporaries thread-local. + */ + padic_t p1; + padic_t p2; + padic_t p3; + /* Used by PackPadic() to build a canonical (reduced) representation. */ + padic_t p4; + /* Scratch rational used by FORM <-> padic/fmpq conversions. */ + fmpq_t q1; +} PADIC_AUX; + +/* + AT.padic_aux_ is allocated by StartPadicSystem() for all threads. + When p-adics are not active this pointer is 0. +*/ +#define PadicAux ((PADIC_AUX *)(AT.padic_aux_)) +#define paux1 (PadicAux->p1) +#define paux2 (PadicAux->p2) +#define paux3 (PadicAux->p3) +#define paux4 (PadicAux->p4) +#define pauxq1 (PadicAux->q1) + +/* + #] Includes : + #[ Helpers : + #[ AllocatePadicPrintBuffer : +*/ +/* + Allocate a conservative print buffer for FLINT's p-adic series printing. + FLINT owns the exact formatting; this just keeps the print buffer large + enough: each coefficient c_n*p^n may need at most + 2*digits_p + digits_exp + 5 + characters as 0<= c_n < p. + We also need room for surrounding parentheses and the string terminator. +*/ +static void AllocatePadicPrintBuffer(LONG ncoeffs, const fmpz_t prime, LONG maxexp) +{ + size_t digits_p = fmpz_sizeinbase(prime,10), digits_exp = 1; + + if ( ncoeffs < 1 ) ncoeffs = 1; + if ( digits_p == 0 ) digits_p = 1; + while ( maxexp >= 10 ) { + maxexp /= 10; + digits_exp++; + } + + if ( AO.padicspace ) M_free(AO.padicspace,"padicspace"); + AO.padicncoeffs = ncoeffs; + AO.padicsize = (LONG)(ncoeffs*(2*digits_p + digits_exp + 5) + 3); + AO.padicspace = (UBYTE *)Malloc1((size_t)AO.padicsize,"padicspace"); +} +/* + #] AllocatePadicPrintBuffer : + #[ InitPadicAux : +*/ +static void InitPadicAux(PADIC_AUX *aux, LONG prec) +{ + padic_init2(aux->p1, (slong)prec); + padic_init2(aux->p2, (slong)prec); + padic_init2(aux->p3, (slong)prec); + padic_init2(aux->p4, (slong)prec); + fmpq_init(aux->q1); +} +/* + #] InitPadicAux : + #[ ClearSinglePadicAux : +*/ +static void ClearSinglePadicAux(PADIC_AUX *aux) +{ + padic_clear(aux->p1); + padic_clear(aux->p2); + padic_clear(aux->p3); + padic_clear(aux->p4); + fmpq_clear(aux->q1); +} +/* + #] ClearSinglePadicAux : + #[ AllocatePadicAux : +*/ +static void AllocatePadicAux(void) +{ +#ifdef WITHPTHREADS + int id, totnum; + totnum = AM.totalnumberofthreads; +#ifdef WITHSORTBOTS + totnum = MaX(2 * AM.totalnumberofthreads - 3, AM.totalnumberofthreads); +#endif + for ( id = 0; id < totnum; id++ ) { + PADIC_AUX *aux; + if ( AB[id]->T.padic_aux_ ) { + ClearSinglePadicAux((PADIC_AUX *)(AB[id]->T.padic_aux_)); + M_free(AB[id]->T.padic_aux_,"AB[id]->T.padic_aux_"); + AB[id]->T.padic_aux_ = 0; + } + aux = Malloc1(sizeof(PADIC_AUX),"AB[id]->T.padic_aux_"); + InitPadicAux(aux,PadicPrecision); + AB[id]->T.padic_aux_ = (void *)aux; + } +#else + PADIC_AUX *aux; + /* + Single-thread build: AT.padic_aux_ is the only instance. + */ + if ( AT.padic_aux_ ) { + ClearSinglePadicAux((PADIC_AUX *)(AT.padic_aux_)); + M_free(AT.padic_aux_,"AT.padic_aux_"); + AT.padic_aux_ = 0; + } + aux = Malloc1(sizeof(PADIC_AUX),"AT.padic_aux_"); + InitPadicAux(aux,PadicPrecision); + AT.padic_aux_ = (void *)aux; +#endif +} +/* + #] AllocatePadicAux : + #[ ClearPadicAux : +*/ +static void ClearPadicAux(void) +{ +#ifdef WITHPTHREADS + int id, totnum; + totnum = AM.totalnumberofthreads; +#ifdef WITHSORTBOTS + totnum = MaX(2 * AM.totalnumberofthreads - 3, AM.totalnumberofthreads); +#endif + for ( id = 0; id < totnum; id++ ) { + if ( AB[id]->T.padic_aux_ ) { + ClearSinglePadicAux((PADIC_AUX *)(AB[id]->T.padic_aux_)); + M_free(AB[id]->T.padic_aux_,"AB[id]->T.padic_aux_"); + AB[id]->T.padic_aux_ = 0; + } + } +#else + if ( AT.padic_aux_ ) { + ClearSinglePadicAux((PADIC_AUX *)(AT.padic_aux_)); + M_free(AT.padic_aux_,"AT.padic_aux_"); + AT.padic_aux_ = 0; + } +#endif +} +/* + #] ClearPadicAux : + #[ StartPadicSystem : + + Initializes (or reinitializes) the single global p-adic context. + Called by #StartPadic from the preprocessor. + + This function: + - clears the previous context if active, + - initializes FLINT's padic_ctx_struct for (p,N), + - allocates per-thread scratch objects (AT.padic_aux_ / AB[id]->T.padic_aux_), + - allocate buffer space for an output string. +*/ +int StartPadicSystem(UBYTE *p, LONG N) +{ + fmpz_t prime; + slong ctxmax; + LONG maxexp; + int error = 0; + fmpz_init(prime); + if ( fmpz_set_str(prime,(char *)p,10) || fmpz_cmp_ui(prime,1) <= 0 ) { + MesPrint("@Illegal prime number in %#StartPadic: %s",p); + error = 1; + } + else if ( fmpz_is_prime(prime) != 1 ) { + MesPrint("@The first parameter in %#StartPadic should be prime: %s",p); + error = 1; + } + if ( error == 0 ) { + if ( PadicActive ) { // Clear the previous padic system + ClearPadicSystem(); + } + PadicPrecision = N; + /* + FLINT contexts cache nonnegative powers p^e, min <= e <= max and + require 0 <= min <= max. For N < 0 an empty cache range is sufficient + and FLINT computes any needed positive powers on demand. + */ + ActivePadicContext = Malloc1(sizeof(padic_ctx_struct),"PadicContext"); + ctxmax = (N > 0) ? (slong)N : 0; + padic_ctx_init(PadicContext,prime,0,ctxmax,PADIC_SERIES); + // Allocate the auxiliary variables + AllocatePadicAux(); + // Allocate buffer space for an output string + maxexp = ABS(N-1); + AllocatePadicPrintBuffer((N > 0) ? N : 1,prime,maxexp); + + PadicActive = 1; + } + fmpz_clear(prime); + return(error); +} +/* + #] StartPadicSystem : + #[ ClearPadicSystem : + + Releases all runtime state associated with p-adic arithmetic, including: + - per-thread aux buffers, + - the FLINT context, + - cached print buffer AO.padicspace. +*/ +void ClearPadicSystem(void) +{ + ClearPadicAux(); + if ( PadicActive ) { + padic_ctx_clear(PadicContext); + } + if ( ActivePadicContext ) { + M_free(ActivePadicContext,"PadicContext"); + ActivePadicContext = 0; + } + PadicActive = 0; + PadicPrecision = 0; + if ( AO.padicspace ) { + M_free(AO.padicspace,"padicspace"); + AO.padicspace = 0; + } + AO.padicsize = 0; + AO.padicncoeffs = 0; +} +/* + #] ClearPadicSystem : + #] Helpers : + #[ Internal p-adic function format : + #[ Explanations : + + p-adic coefficients are stored in a special function padic_ whose + arguments contain the relevant p-adic information from the FLINT library: + padic_(v, N, u) + where + - v is the p-adic valuation, + - N is the precision stored with this coefficient, + - u is the FLINT unit. + + The prime p is not stored in the function. It is supplied by the active + p-adic context configured through #StartPadic / #EndPadic. The unit + u is stored in the numerator of a normal FORM integer with denominator 1. + Note: we could store p in the denominator as p cannot divide u. + + #] Explanations : + #[ TestPadic : + + Checks whether fun has the shape of a legal padic_ with its arguments: + padic_(v,N,u) + + Return value: + - 1 if fun is a syntactically valid padic_ function, + - 0 otherwise. +*/ +int TestPadic(WORD *fun) +{ + WORD *f, *fstop; + WORD nargs, nnum, i; + + f = fun + FUNHEAD; + fstop = fun + fun[1]; + nargs = 0; + while ( f < fstop ) { + nargs++; + NEXTARG(f); + } + if ( nargs != 3 || f != fstop ) return(0); + + f = fun + FUNHEAD; + /* v: signed LONG argument */ + if ( *f == -SNUMBER ) { + f += 2; + } + else { + if ( *f != ARGHEAD+6 ) return(0); + if ( ABS(f[ARGHEAD+5]) != 5 ) return(0); + if ( f[ARGHEAD+3] != 1 ) return(0); + if ( f[ARGHEAD+4] != 0 ) return(0); + f += *f; + } + /* N: signed LONG argument. FLINT supports absolute precision at, + above, and below O(p^0). */ + if ( *f == -SNUMBER ) { + f += 2; + } + else { + if ( *f != ARGHEAD+6 ) return(0); + if ( ABS(f[ARGHEAD+5]) != 5 ) return(0); + if ( f[ARGHEAD+3] != 1 ) return(0); + if ( f[ARGHEAD+4] != 0 ) return(0); + f += *f; + } + /* u: FORM integer (-SNUMBER or long integer n/1). */ + if ( *f == -SNUMBER ) { + f += 2; + } + else { + nnum = (WORD)((ABS(f[*f-1])-1)/2); + if ( nnum <= 0 ) return(0); + if ( *f != ARGHEAD + 2*nnum + 2 ) return(0); + if ( f[ARGHEAD] != 2*nnum + 2 ) return(0); + /* Denominator must be exactly 1. */ + f += ARGHEAD + nnum + 1; + if ( f[0] != 1 ) return(0); + for ( i = 1; i < nnum; i++ ) { + if ( f[i] != 0 ) return(0); + } + } + return(1); +} +/* + #] TestPadic : + #[ UnpackPadic : + + Converts the arguments of a padic_ function back into FLINT's padic_t. + We shouldn't come here if the p-adic system is turned off. + + Return value: + - 0 on success, + - -1 if fun is not a valid internal p-adic function. +*/ +static int UnpackPadic(padic_t out, WORD *fun) +{ + WORD *f; + ULONG x; + + if ( !PadicActive ) { + MLOCK(ErrorMessageLock); + MesPrint("Illegal attempt at using a padic_ function without proper startup."); + MesPrint("Please use %#StartPadic

, first."); + MUNLOCK(ErrorMessageLock); + Terminate(-1); + } + + if ( TestPadic(fun) == 0 ) return(-1); + + /* + Read the argument triplet in order: v, N, u. + */ + f = fun + FUNHEAD; + + /* + TestPadic() already validated the arguments of padic_, + so we can safely decode v, N and u here. + */ + if ( *f == -SNUMBER ) { + padic_val(out) = (slong)f[1]; + f += 2; + } + else { + x = ((ULONG)(UWORD)f[ARGHEAD+2] << BITSINWORD) + (UWORD)f[ARGHEAD+1]; + padic_val(out) = (f[ARGHEAD+5] < 0) ? -(slong)x : (slong)x; + f += *f; + } + + if ( *f == -SNUMBER ) { + padic_prec(out) = (slong)f[1]; + f += 2; + } + else { + x = ((ULONG)(UWORD)f[ARGHEAD+2] << BITSINWORD) + (UWORD)f[ARGHEAD+1]; + padic_prec(out) = (f[ARGHEAD+5] < 0) ? -(slong)x : (slong)x; + f += *f; + } + + /* + Decode the unit u. + */ + if ( *f == -SNUMBER ) { + fmpz_set_si(padic_unit(out),(slong)(f[1])); + } + else { + WORD size, nnum; + size = f[*f-1]; + nnum = (WORD)((ABS(size)-1)/2); + flint_fmpz_set_form(padic_unit(out),(UWORD *)(f+ARGHEAD+1), + (size < 0) ? -nnum : nnum); + } + + padic_reduce(out,PadicContext); + return(0); +} +/* + #] UnpackPadic : + #[ PackPadic : + + Converts a FLINT's padic_t into FORM's internal padic notation: + padic_(v,N,u) + + Return value: + - the number of WORDs written to fun. +*/ +static int PackPadic(WORD *fun, padic_t in) +{ + WORD *t; + LONG v, N, small; + ULONG x; + WORD *arg; + WORD nnum, nabs, i; + GETIDENTITY + + /* + Normalize first so the (v,N,u) triplet is canonical. + */ + padic_set(paux4,in,PadicContext); + padic_reduce(paux4,PadicContext); + v = (LONG)padic_val(paux4); + N = (LONG)padic_prec(paux4); + + /* + We now fill the function with the three arguments: + valuation v, precision N, and unit u. + */ + t = fun; + *t++ = PADICFUN; + t++; // fun[1] (function length) is filled at the end. + FILLFUN(t); + + /* + Pack valuation v: + compact -SNUMBER when it fits in WORD, otherwise as a two-word + long integer argument. + */ + if ( v >= WORD_MIN_VALUE && v <= WORD_MAX_VALUE ) { + *t++ = -SNUMBER; + *t++ = (WORD)v; + } + else { + x = ( v < 0 ) ? (ULONG)(-v) : (ULONG)v; + *t++ = ARGHEAD + 6; + *t++ = 0; + FILLARG(t); + *t++ = 6; + *t++ = (UWORD)x; + *t++ = (UWORD)(x >> BITSINWORD); + *t++ = 1; + *t++ = 0; + *t++ = (v < 0) ? -5 : 5; + } + /* + Pack precision N. + */ + if ( N >= WORD_MIN_VALUE && N <= WORD_MAX_VALUE ) { + *t++ = -SNUMBER; + *t++ = (WORD)N; + } + else { + x = ( N < 0 ) ? (ULONG)(-N) : (ULONG)N; + *t++ = ARGHEAD + 6; + *t++ = 0; + FILLARG(t); + *t++ = 6; + *t++ = (UWORD)x; + *t++ = (UWORD)(x >> BITSINWORD); + *t++ = 1; + *t++ = 0; + *t++ = (N < 0) ? -5 : 5; + } + /* + Pack unit u in the numerator of a FORM integer. + use -SNUMBER for small values, otherwise FORM long-integer. + */ + if ( fmpz_fits_si(padic_unit(paux4)) + && (small = (LONG)fmpz_get_si(padic_unit(paux4)), + small >= WORD_MIN_VALUE && small <= WORD_MAX_VALUE) ) { + *t++ = -SNUMBER; + *t++ = (WORD)small; + } + else { + arg = t + ARGHEAD + 1; + // Numerator + nnum = flint_fmpz_get_form(padic_unit(paux4),arg); + nabs = ABS(nnum); + // Denominator + arg += nabs; + *arg++ = 1; + for ( i = 1; i < nabs; i++ ) *arg++ = 0; + *arg++ = ( nnum < 0 ) ? -(2*nabs+1) : (2*nabs+1); + // Arghead etc. + *t++ = ARGHEAD + 2*nabs + 2; + *t++ = 0; + FILLARG(t); + *t++ = 2*nabs + 2; + t = arg; + } + fun[1] = t - fun; + return(fun[1]); +} +/* + #] PackPadic : + #] Internal p-adic function format : + #[ Rekenen : + #[ FormRatToFmpq : + + Converts the internal FORM rational coefficient encoding (formrat/ratsize) + to a FLINT rational fmpq. + + Maybe this should be moved to flintinterface.cc? +*/ +static void FormRatToFmpq(fmpq_t result, UWORD *formrat, WORD ratsize) +{ + WORD nnum, nden; + UWORD *num, *den; + int sign = 1; + + if ( ratsize < 0 ) { + sign = -1; + ratsize = -ratsize; + } + + nnum = nden = (ratsize-1)/2; + num = formrat; + den = formrat + nnum; + /* Remove padding */ + while ( nnum > 0 && num[nnum-1] == 0 ) nnum--; + while ( nden > 0 && den[nden-1] == 0 ) nden--; + + if ( nnum > 0 ) { + flint_fmpz_set_form(fmpq_numref(result),num,sign*nnum); + } + else { + fmpz_zero(fmpq_numref(result)); + } + flint_fmpz_set_form(fmpq_denref(result),den,nden); + + fmpq_canonicalise(result); +} +/* + #] FormRatToFmpq : + #[ MulRatToPadic : + + Multiplies an existing `padic_` coefficient by a FORM rational coefficient. + This is used by Normalize() when a term contains both a numeric coefficient + and a `padic_` function. + + Return value: + - 0 on success, with a non-zero product packed into outfun, + - 1 if the product is p-adic zero, + - -1 on runtime/validation failure. +*/ +int MulRatToPadic(PHEAD WORD *outfun, WORD *infun, UWORD *formrat, WORD nrat) +{ + if ( !PadicActive ) return(-1); + if ( UnpackPadic(paux1,infun) ) return(-1); + FormRatToFmpq(pauxq1,formrat,nrat); + padic_set_fmpq(paux2,pauxq1,PadicContext); + padic_mul(paux3,paux1,paux2,PadicContext); + PackPadic(outfun,paux3); + if ( padic_is_zero(paux3) ) return(1); + return(0); +} +/* + #] MulRatToPadic : + #[ MulPadics : + + Multiplies two padic_ functions (fun1 * fun2) and stores the + product as a new padic_ function in fun3. + + Return value: + - 0 on success, with a non-zero product packed into fun3, + - 1 if the product is p-adic zero, + - -1 on runtime/validation failure. +*/ +int MulPadics(PHEAD WORD *fun3, WORD *fun1, WORD *fun2) +{ + if ( !PadicActive ) return(-1); + if ( UnpackPadic(paux1,fun1) ) return(-1); + if ( UnpackPadic(paux2,fun2) ) return(-1); + padic_mul(paux3,paux1,paux2,PadicContext); + PackPadic(fun3,paux3); + if ( padic_is_zero(paux3) ) return(1); + return(0); +} +/* + #] MulPadics : + #[ DivPadics : + + Divides two padic_ functions (fun1 / fun2) and stores the + quotient as a new padic_ function in fun3. + Division by zero is a fatal runtime error. + + Return value: + - 0 on success, with a non-zero quotient packed into fun3, + - 1 if the quotient is p-adic zero, + - -1 on runtime/validation failure. +*/ +int DivPadics(PHEAD WORD *fun3, WORD *fun1, WORD *fun2) +{ + if ( !PadicActive ) return(-1); + if ( UnpackPadic(paux1,fun1) ) return(-1); + if ( UnpackPadic(paux2,fun2) ) return(-1); + if ( padic_is_zero(paux2) ) { + MLOCK(ErrorMessageLock); + MesPrint("Division by zero in p-adic arithmetic."); + MUNLOCK(ErrorMessageLock); + Terminate(-1); + return(-1); + } + padic_div(paux3,paux1,paux2,PadicContext); + PackPadic(fun3,paux3); + if ( padic_is_zero(paux3) ) return(1); + return(0); +} +/* + #] DivPadics : + #[ InvPadic : + + Computes the inverse 1/fun and stores the result as a new padic_ function + in outfun. Division by zero is a fatal runtime error. + + Return value: + - 0 on success, with a non-zero inverse packed into outfun, + - 1 if the inverse is p-adic zero, + - -1 on runtime/validation failure. +*/ +int InvPadic(PHEAD WORD *outfun, WORD *fun) +{ + if ( !PadicActive ) return(-1); + if ( UnpackPadic(paux1,fun) ) return(-1); + if ( padic_is_zero(paux1) ) { + MLOCK(ErrorMessageLock); + MesPrint("Division by zero in p-adic arithmetic."); + MUNLOCK(ErrorMessageLock); + Terminate(-1); + return(-1); + } + padic_inv(paux3,paux1,PadicContext); + PackPadic(outfun,paux3); + if ( padic_is_zero(paux3) ) return(1); + return(0); +} +/* + #] InvPadic : + #[ PadicReconstruct : + + Reconstructs a FORM rational coefficient from a p-adic coefficient. + + This mirrors the FORM-level recipe + + FromPadic, padic; + id padic(v?,N?,u?) = makerational_(u,`$prime'^(N-v))*$prime^v; + + but calls the same low-level reconstruction routines used by makerational_ + directly. The output is a full FORM coefficient, including the final signed + length word. A return value of -1 is non-fatal: no unique reconstruction was + found or an intermediate value does not fit local number storage, so + PadicToRat should leave the original padic_ untouched. +*/ +static int PadicReconstruct(PHEAD UWORD *out, WORD *nratout, WORD *fun) +{ + fmpz_t modulus, ppower; + UWORD *u, *mod, *formpower; + WORD *arg; + WORD nu, nmod, npower, nred, i; + ULONG x; + slong v, N, modexp; + int retval = -1; + + if ( TestPadic(fun) == 0 ) return(-1); + + // Read the argument triplet in order: v, N, u. + arg = fun + FUNHEAD; + if ( *arg == -SNUMBER ) { + v = (slong)arg[1]; + arg += 2; + } + else { + x = ((ULONG)(UWORD)arg[ARGHEAD+2] << BITSINWORD) + (UWORD)arg[ARGHEAD+1]; + v = (arg[ARGHEAD+5] < 0) ? -(slong)x : (slong)x; + arg += *arg; + } + if ( *arg == -SNUMBER ) { + N = (slong)arg[1]; + arg += 2; + } + else { + x = ((ULONG)(UWORD)arg[ARGHEAD+2] << BITSINWORD) + (UWORD)arg[ARGHEAD+1]; + N = (arg[ARGHEAD+5] < 0) ? -(slong)x : (slong)x; + arg += *arg; + } + + if ( *arg == -SNUMBER && arg[1] == 0 ) { + out[0] = 0; + out[1] = 1; + *nratout = INCLENG(1); + out[2] = (UWORD)ABS(*nratout); + return(0); + } + modexp = N - v; + if ( modexp <= 0 ) return(-1); + + fmpz_init(modulus); + fmpz_init(ppower); + u = NumberMalloc("PadicReconstruct"); + mod = NumberMalloc("PadicReconstruct"); + formpower = NumberMalloc("PadicReconstruct"); + + if ( *arg == -SNUMBER ) { + nu = (arg[1] < 0) ? -1 : 1; + u[0] = (arg[1] < 0) ? (UWORD)(-arg[1]) : (UWORD)arg[1]; + } + else { + nu = (WORD)((ABS(arg[*arg-1])-1)/2); + if ( arg[*arg-1] < 0 ) nu = -nu; + for ( i = 0; i < ABS(nu); i++ ) u[i] = (UWORD)arg[ARGHEAD+1+i]; + } + + /* mod = modulus = p^(N-v). */ + fmpz_pow_ui(modulus,PadicContext->p,(ulong)modexp); + nmod = flint_fmpz_get_form(modulus,(WORD *)mod); + + /* + Match makerational_: rational reconstruction works with the residue + class u mod p^(N-v), so reduce u in place when it lies outside the + modulus range. The quotient is not needed. + */ + if ( BigLong(u,ABS(nu),mod,nmod) >= 0 ) { + UWORD *quotient = NumberMalloc("PadicReconstruct"); + UWORD *remainder = NumberMalloc("PadicReconstruct"); + WORD nquot, nrem; + + if ( DivLong(u,nu,mod,nmod,quotient,&nquot,remainder,&nrem) ) { + NumberFree(remainder,"PadicReconstruct"); + NumberFree(quotient,"PadicReconstruct"); + goto ClearAndReturn; + } + for ( i = 0; i < ABS(nrem); i++ ) u[i] = remainder[i]; + nu = nrem; + NumberFree(quotient,"PadicReconstruct"); + NumberFree(remainder,"PadicReconstruct"); + if ( nu == 0 ) { + out[0] = 0; + out[1] = 1; + nred = 1; + goto StoreCoefficient; + } + } + if ( ABS(nu) == 1 && nmod == 1 + && u[0] <= (UWORD)WORD_MAX_VALUE + && mod[0] <= (UWORD)WORD_MAX_VALUE ) { + WORD num, den; + WORD sign = (nu < 0) ? -1 : 1; + + if ( MakeRational((WORD)u[0],(WORD)mod[0],&num,&den) ) goto ClearAndReturn; + if ( sign < 0 ) num = -num; + if ( num < 0 ) { + out[0] = (UWORD)(-num); + nred = -1; + } + else { + out[0] = (UWORD)num; + nred = 1; + } + out[1] = (UWORD)den; + } + else { + if ( MakeLongRational(BHEAD u,nu,mod,nmod,out,&nred) ) { + goto ClearAndReturn; + } + } + + /* + The reconstruction above returns only the rational unit. Now apply the + p-adic valuation with FORM's rational arithmetic so the coefficient is + normalized in the same way as makerational_ output. + */ + if ( v > 0 ) { + fmpz_pow_ui(ppower,PadicContext->p,(ulong)v); + npower = flint_fmpz_get_form(ppower,(WORD *)formpower); + if ( Mully(BHEAD out,&nred,formpower,npower) ) goto ClearAndReturn; + } + else if ( v < 0 ) { + fmpz_pow_ui(ppower,PadicContext->p,(ulong)(-v)); + npower = flint_fmpz_get_form(ppower,(WORD *)formpower); + if ( Divvy(BHEAD out,&nred,formpower,npower) ) goto ClearAndReturn; + } + +StoreCoefficient: + *nratout = INCLENG(nred); + out[2*ABS(nred)] = (UWORD)ABS(*nratout); + retval = 0; + +ClearAndReturn: + NumberFree(formpower,"PadicReconstruct"); + NumberFree(mod,"PadicReconstruct"); + NumberFree(u,"PadicReconstruct"); + fmpz_clear(ppower); + fmpz_clear(modulus); + return(retval); +} +/* + #] PadicReconstruct : + #] Rekenen : + #[ Printing : + #[ PrintPadic : + + Formats a padic_ function for printing in series mode as: + (c_v*p^v + ... + c_{N-1}*p^(N-1)) + with v the valuation and N the precision. + + The resulting C string is stored in AO.padicspace and the return value + is the string length, or 0 if no printable p-adic string was produced. +*/ +int PrintPadic(WORD *fun) +{ + GETIDENTITY + char *out, *series; + slong v, N; + size_t n; + LONG ncoeffs, maxexp; + + if ( !PadicActive ) return(0); + if ( UnpackPadic(paux1,fun) ) return(0); + + v = padic_val(paux1); + N = padic_prec(paux1); + + if ( N <= v ) ncoeffs = 1; + else ncoeffs = (LONG)(N-v); + + /* Grow the print buffer when it needs to accomodate more coefficients. */ + if ( ncoeffs > AO.padicncoeffs ) { + maxexp = MaX(ABS((LONG)v),(LONG)N); + AllocatePadicPrintBuffer(ncoeffs,PadicContext->p,maxexp); + } + + out = (char *)AO.padicspace; + out[0] = '('; + series = padic_get_str(out+1,paux1,PadicContext); + if ( series == 0 ) return(0); + n = strlen(series); + + out[1+n] = ')'; + out[2+n] = 0; + return((int)(n+2)); +} +/* + #] PrintPadic : + #] Printing : + #[ Compiler/runtime statements : + #[ CoToPadic : + + Compiler front-end for the `ToPadic;` statement. + This only validates syntax and records the action; the actual conversion is + performed on terms at execution time by ToPadic(). +*/ +int CoToPadic(UBYTE *s) +{ + if ( !PadicActive ) { + MesPrint("&Illegal attempt to convert to padic_ without activating p-adic numbers."); + MesPrint("&Forgotten %#startpadic instruction?"); + return(1); + } + while ( *s == ' ' || *s == ',' || *s == '\t' ) s++; + if ( *s ) { + MesPrint("&Illegal argument(s) in ToPadic statement: '%s'",s); + return(1); + } + Add2Com(TYPETOPADIC); + return(0); +} +/* + #] CoToPadic : + #[ CoPadicToRat : + + Compiler front-end for the `PadicToRat;` statement. + The runtime action is implemented by PadicToRat(). +*/ +int CoPadicToRat(UBYTE *s) +{ + if ( !PadicActive ) { + MesPrint("&Illegal attempt to convert from padic_ without activating p-adic numbers."); + MesPrint("&Forgotten %#startpadic instruction?"); + return(1); + } + while ( *s == ' ' || *s == ',' || *s == '\t' ) s++; + if ( *s ) { + MesPrint("&Illegal argument(s) in PadicToRat statement: '%s'",s); + return(1); + } + Add2Com(TYPEPADICTORAT); + return(0); +} +/* + #] CoPadicToRat : + #[ CoFromPadic : + + Compiler front-end for the `FromPadic, f;` statement. + The target must be a declared regular function. +*/ +int CoFromPadic(UBYTE *s) +{ + WORD numfun; + int type; + UBYTE *t, c; + + while ( *s == ' ' || *s == ',' || *s == '\t' ) s++; + t = SkipAName(s); + if ( t == 0 ) goto syntaxerror; + c = *t; *t = 0; + type = GetName(AC.varnames,s,&numfun,NOAUTO); + *t = c; + if ( type != CFUNCTION || functions[numfun].spec != 0 + || numfun + FUNCTION == PADICFUN ) goto syntaxerror; + while ( *t == ' ' || *t == ',' || *t == '\t' ) t++; + if ( *t ) goto syntaxerror; + Add3Com(TYPEFROMPADIC,numfun+FUNCTION); + return(0); + +syntaxerror: + MesPrint("&FromPadic statement needs one declared regular function for its argument"); + return(1); +} +/* + #] CoFromPadic : + #[ ToPadic : + + Runtime implementation of `ToPadic;`. + + This replaces the current rational coefficient by an explicit padic_ function + and sets the rational coefficient to 1/1. The sign is absorbed into the + p-adic value. If the coefficient is already 1/1 and the term already ends + with a proper padic_ function we are done immediately. +*/ +int ToPadic(PHEAD WORD *term, WORD level) +{ + GETBIDENTITY + WORD *t, *tstop, nsize, ncoef; + + if ( !PadicActive ) return(1); + + t = term + *term; + ncoef = t[-1]; + nsize = ABS(ncoef); + tstop = t - nsize; + + /* If there is already a proper padic_, with rational coefficient 1/1 we are done. */ + if ( ncoef == 3 && t[-2] == 1 && t[-3] == 1 ) { + t = term + 1; + while ( t < tstop ) { + if ( *t == PADICFUN && (t+t[1] == tstop) && TestPadic(t) ) { + return(Generator(BHEAD term,level)); + } + t += t[1]; /* advance to the next function in the term */ + } + } + + FormRatToFmpq(pauxq1,(UWORD *)tstop,ncoef); + padic_set_fmpq(paux1,pauxq1,PadicContext); + if ( padic_is_zero(paux1) ) return(0); + // Overwrite the old rational coefficient of the term by padic_(v,N,u) + // and append rational coefficient 1/1. + PackPadic(tstop,paux1); + tstop += tstop[1]; /* advance past the newly written padic_ */ + *tstop++ = 1; + *tstop++ = 1; + *tstop++ = 3; + *term = tstop - term; + AT.WorkPointer = tstop; + return(Generator(BHEAD term,level)); +} +/* + #] ToPadic : + #[ PadicToRat : + + Runtime implementation of `PadicToRat;`. + + This finds a terminal `padic_` coefficient record and converts it back to a + FORM rational coefficient by rational reconstruction. If no reconstruction + exists at the current precision, the padic_ coefficient is left untouched. +*/ +int PadicToRat(PHEAD WORD *term, WORD level) +{ + GETBIDENTITY + static int warnflag = 1; + WORD *tstop, *t, *stop, *from, nsize, nsign, ncoef, i; + UWORD *rat; + + if ( !PadicActive ) return(1); + + tstop = term + *term; + nsize = ABS(tstop[-1]); + nsign = tstop[-1] < 0 ? -1 : 1; + tstop -= nsize; + t = term + 1; + /* + The term must end in a single proper padic_ function, followed by a unit + coefficient 1/1 (the sign is carried separately in tstop[-1]). + */ + while ( t < tstop ) { + if ( *t == PADICFUN && t + t[1] == tstop && TestPadic(t) + && nsize == 3 && tstop[0] == 1 && tstop[1] == 1 ) break; + t += t[1]; + } + if ( t < tstop ) { + rat = (UWORD *)TermMalloc("PadicToRat"); + if ( PadicReconstruct(BHEAD rat,&ncoef,t) ) { + if ( warnflag ) { + MLOCK(ErrorMessageLock); + if ( warnflag ) { + MesPrint("%w Warning: p-adic coefficient could not be reconstructed in PadicToRat"); + warnflag = 0; + } + MUNLOCK(ErrorMessageLock); + } + TermFree((WORD *)rat,"PadicToRat"); + return(Generator(BHEAD term,level)); + } + if ( rat[0] == 0 && rat[1] == 1 && ncoef == 3 ) { + TermFree((WORD *)rat,"PadicToRat"); + return(0); + } + stop = (WORD *)(((UBYTE *)term) + AM.MaxTer); + if ( t + ABS(ncoef) > stop ) { + MLOCK(ErrorMessageLock); + MesPrint("Term too complex after p-adic to rational conversion. MaxTermSize = %10l", + AM.MaxTer/sizeof(WORD)); + MUNLOCK(ErrorMessageLock); + TermFree((WORD *)rat,"PadicToRat"); + Terminate(-1); + return(1); + } + from = (WORD *)rat; + i = ABS(ncoef); + NCOPY(t,from,i) + t[-1] = ncoef*nsign; + *term = t - term; + TermFree((WORD *)rat,"PadicToRat"); + } + return(Generator(BHEAD term,level)); +} +/* + #] PadicToRat : + #[ FromPadic : + + Runtime implementation of `FromPadic, f;`. + Each syntactically valid padic_(v,N,u) at the current execution level is + rewritten in place to f(v,N,u). +*/ +int FromPadic(PHEAD WORD *term, WORD level, WORD target) +{ + GETBIDENTITY + WORD *t, *tstop; + + tstop = term + *term; + tstop -= ABS(tstop[-1]); + t = term + 1; + while ( t < tstop ) { + if ( *t == PADICFUN && TestPadic(t) ) *t = target; + t += t[1]; + } + return(Generator(BHEAD term,level)); +} +/* + #] FromPadic : + #] Compiler/runtime statements : + #[ Sorting : + #[ AddWithPadic : + + Sort helper used when two otherwise identical terms differ only in their + coefficient and that coefficient involves padic_. + + Compare1() in sort.c sets AT.SortPadicMode to indicate which term(s) carry a + padic_ coefficient record: + - 3: both terms end with a proper padic_ record + - 1: only *ps1 has padic_, *ps2 has a rational coefficient + - 2: only *ps2 has padic_, *ps1 has a rational coefficient + + AddWithPadic computes the p-adic sum and rewrites *ps1 to contain exactly one + padic_ record and a unit coefficient 1/1. It tries to reuse term space in + place and otherwise allocates in the sort buffer. +*/ +int AddWithPadic(PHEAD WORD **ps1, WORD **ps2) +{ + GETBIDENTITY + SORTING *S = AT.SS; + WORD *coef1, *coef2, size1, size2, *fun1, *fun2, *fun3; + WORD *s1, *s2, *t1, *t2, i, j, jj; + + if ( !PadicActive ) return(0); + + s1 = *ps1; + s2 = *ps2; + coef1 = s1 + *s1; size1 = coef1[-1]; coef1 -= ABS(size1); + coef2 = s2 + *s2; size2 = coef2[-1]; coef2 -= ABS(size2); + + if ( AT.SortPadicMode == 3 ) { + /* Both coefficients are padic_: unpack and apply external +/- sign. */ + fun1 = s1+1; while ( fun1 < coef1 && fun1[0] != PADICFUN ) fun1 += fun1[1]; + fun2 = s2+1; while ( fun2 < coef2 && fun2[0] != PADICFUN ) fun2 += fun2[1]; + UnpackPadic(paux1,fun1); + if ( size1 < 0 ) padic_neg(paux1,paux1,PadicContext); + UnpackPadic(paux2,fun2); + if ( size2 < 0 ) padic_neg(paux2,paux2,PadicContext); + } + else if ( AT.SortPadicMode == 1 ) { + /* First coefficient is padic_, second is rational. */ + fun1 = s1+1; while ( fun1 < coef1 && fun1[0] != PADICFUN ) fun1 += fun1[1]; + UnpackPadic(paux1,fun1); + if ( size1 < 0 ) padic_neg(paux1,paux1,PadicContext); + fun2 = coef2; + FormRatToFmpq(pauxq1,(UWORD *)coef2,size2); + padic_set_fmpq(paux2,pauxq1,PadicContext); + } + else if ( AT.SortPadicMode == 2 ) { + /* Second coefficient is padic_, first is rational. */ + fun2 = s2+1; while ( fun2 < coef2 && fun2[0] != PADICFUN ) fun2 += fun2[1]; + UnpackPadic(paux2,fun2); + if ( size2 < 0 ) padic_neg(paux2,paux2,PadicContext); + fun1 = coef1; + FormRatToFmpq(pauxq1,(UWORD *)coef1,size1); + padic_set_fmpq(paux1,pauxq1,PadicContext); + } + else { + MLOCK(ErrorMessageLock); + MesPrint("Illegal value %d for AT.SortPadicMode in AddWithPadic.",AT.SortPadicMode); + MUNLOCK(ErrorMessageLock); + Terminate(-1); + return(0); + } + + padic_add(paux3,paux1,paux2,PadicContext); + if ( padic_is_zero(paux3) ) { + /* Terms cancel. */ + *ps1 = *ps2 = 0; + AT.SortPadicMode = 0; + return(0); + } + + fun3 = TermMalloc("AddWithPadic"); + PackPadic(fun3,paux3); + + if ( AT.SortPadicMode == 3 ) { + /* Prefer overwriting an existing padic_ record in-place if it fits. */ + if ( fun1[1] == fun3[1] ) { +Over1: + i = fun3[1]; t1 = fun1; t2 = fun3; NCOPY(t1,t2,i); + *t1++ = 1; *t1++ = 1; *t1++ = 3; + *s1 = t1-s1; goto Finished; + } + else if ( fun2[1] == fun3[1] ) { +Over2: + i = fun3[1]; t1 = fun2; t2 = fun3; NCOPY(t1,t2,i); + *t1++ = 1; *t1++ = 1; *t1++ = 3; + *s2 = t1-s2; *ps1 = s2; goto Finished; + } + else if ( fun1[1] >= fun3[1] ) goto Over1; + else if ( fun2[1] >= fun3[1] ) goto Over2; + } + else if ( AT.SortPadicMode == 1 ) { + if ( fun1[1] >= fun3[1] ) goto Over1; + else if ( fun3[1]+3 <= ABS(size2) ) goto Over2; + } + else if ( AT.SortPadicMode == 2 ) { + if ( fun2[1] >= fun3[1] ) goto Over2; + else if ( fun3[1]+3 <= ABS(size1) ) goto Over1; + } + + jj = fun1-s1; + /* Required space: prefix + padic_ record + coefficient words (1,1,3). */ + j = jj+fun3[1]+3; + if ( (S->sFill + j) >= S->sTop2 ) { + /* Make space in the sort buffer and refresh pointers. */ + GarbHand(); + s1 = *ps1; + fun1 = s1+jj; + } + t1 = S->sFill; + for ( i = 0; i < jj; i++ ) *t1++ = s1[i]; + i = fun3[1]; s1 = fun3; NCOPY(t1,s1,i); + *t1++ = 1; *t1++ = 1; *t1++ = 3; + *ps1 = S->sFill; + **ps1 = t1-*ps1; + S->sFill = t1; + +Finished: + *ps2 = 0; + TermFree(fun3,"AddWithPadic"); + AT.SortPadicMode = 0; + if ( **ps1 > AM.MaxTer/((LONG)(sizeof(WORD))) ) { + MLOCK(ErrorMessageLock); + MesPrint("Term too complex after p-adic addition in sort. MaxTermSize = %10l", + AM.MaxTer/sizeof(WORD)); + MUNLOCK(ErrorMessageLock); + Terminate(-1); + } + return(1); +} +/* + #] AddWithPadic : + #[ MergeWithPadic : + + Variant of AddWithPadic used during patch merging: it computes the p-adic sum + and rewrites term1 in-place if possible, otherwise it shifts/overwrites + memory so that *interm1 points to the updated term. +*/ +int MergeWithPadic(PHEAD WORD **interm1, WORD **interm2) +{ + GETBIDENTITY + WORD *coef1, *coef2, size1, size2, *fun1, *fun2, *fun3, *tt; + WORD jj, *t1, *t2, i, *term1 = *interm1, *term2 = *interm2; + int retval = 0; + + if ( !PadicActive ) return(0); + + coef1 = term1+*term1; size1 = coef1[-1]; coef1 -= ABS(size1); + coef2 = term2+*term2; size2 = coef2[-1]; coef2 -= ABS(size2); + if ( AT.SortPadicMode == 3 ) { + fun1 = term1+1; while ( fun1 < coef1 && fun1[0] != PADICFUN ) fun1 += fun1[1]; + fun2 = term2+1; while ( fun2 < coef2 && fun2[0] != PADICFUN ) fun2 += fun2[1]; + UnpackPadic(paux1,fun1); + if ( size1 < 0 ) padic_neg(paux1,paux1,PadicContext); + UnpackPadic(paux2,fun2); + if ( size2 < 0 ) padic_neg(paux2,paux2,PadicContext); + } + else if ( AT.SortPadicMode == 1 ) { + fun1 = term1+1; while ( fun1 < coef1 && fun1[0] != PADICFUN ) fun1 += fun1[1]; + UnpackPadic(paux1,fun1); + if ( size1 < 0 ) padic_neg(paux1,paux1,PadicContext); + fun2 = coef2; + FormRatToFmpq(pauxq1,(UWORD *)coef2,size2); + padic_set_fmpq(paux2,pauxq1,PadicContext); + } + else if ( AT.SortPadicMode == 2 ) { + fun2 = term2+1; while ( fun2 < coef2 && fun2[0] != PADICFUN ) fun2 += fun2[1]; + fun1 = coef1; + FormRatToFmpq(pauxq1,(UWORD *)coef1,size1); + padic_set_fmpq(paux1,pauxq1,PadicContext); + UnpackPadic(paux2,fun2); + if ( size2 < 0 ) padic_neg(paux2,paux2,PadicContext); + } + else { + MLOCK(ErrorMessageLock); + MesPrint("Illegal value %d for AT.SortPadicMode in MergeWithPadic.",AT.SortPadicMode); + MUNLOCK(ErrorMessageLock); + Terminate(-1); + return(0); + } + + padic_add(paux3,paux1,paux2,PadicContext); + if ( padic_is_zero(paux3) ) { + AT.SortPadicMode = 0; + return(0); + } + + fun3 = TermMalloc("MergeWithPadic"); + PackPadic(fun3,paux3); + if ( AT.SortPadicMode == 3 ) { + if ( fun1[1] + ABS(size1) == fun3[1] + 3 ) { +OnTopOf1: + /* The new (padic_ + 1/1) fits exactly on top of the old suffix. */ + t1 = fun3; t2 = fun1; + for ( i = 0; i < fun3[1]; i++ ) *t2++ = *t1++; + *t2++ = 1; *t2++ = 1; *t2++ = 3; + retval = 1; + } + else if ( fun1[1] + ABS(size1) > fun3[1] + 3 ) { +Shift1: + /* There is slack in term1; shift the tail down and rewrite in place. */ + t2 = term1 + *term1; tt = t2; + *--t2 = 3; *--t2 = 1; *--t2 = 1; + t1 = fun3 + fun3[1]; + for ( i = 0; i < fun3[1]; i++ ) *--t2 = *--t1; + t1 = fun1; + while ( t1 > term1 ) *--t2 = *--t1; + *t2 = tt-t2; term1 = t2; + retval = 1; + } + else { + jj = fun3[1]-fun1[1]+3-ABS(size1); +Over1: + /* term1 needs to grow: move the start pointer back by jj words. */ + t2 = term1-jj; t1 = term1; + while ( t1 < fun1 ) *t2++ = *t1++; + term1 -= jj; + *term1 += jj; + for ( i = 0; i < fun3[1]; i++ ) *t2++ = fun3[i]; + *t2++ = 1; *t2++ = 1; *t2++ = 3; + retval = 1; + } + } + else if ( AT.SortPadicMode == 1 ) { + if ( fun1[1] + ABS(size1) == fun3[1] + 3 ) goto OnTopOf1; + else if ( fun1[1] + ABS(size1) > fun3[1] + 3 ) goto Shift1; + else { + jj = fun3[1]-fun1[1]+3-ABS(size1); + goto Over1; + } + } + else { + if ( fun3[1] + 3 == ABS(size1) ) goto OnTopOf1; + else if ( fun3[1] + 3 < ABS(size1) ) goto Shift1; + else { + jj = fun3[1]+3-ABS(size1); + goto Over1; + } + } + *interm1 = term1; + TermFree(fun3,"MergeWithPadic"); + AT.SortPadicMode = 0; + return(retval); +} +/* + #] MergeWithPadic : + #] Sorting : +*/ diff --git a/sources/pre.c b/sources/pre.c index cee07818..ca6c764f 100644 --- a/sources/pre.c +++ b/sources/pre.c @@ -76,6 +76,9 @@ static KEYWORD precommands[] = { ,{"endif" , DoEndif , 0, 0} ,{"endinside" , DoEndInside , 0, 0} ,{"endnamespace" , DoEndNamespace , 0, 0} +#ifdef WITHFLINT + ,{"endpadic" , DoEndPadic , 0, 0} +#endif ,{"endprocedure" , DoEndprocedure , 0, 0} ,{"endswitch" , DoPreEndSwitch , 0, 0} ,{"exchange" , DoPreExchange , 0, 0} @@ -113,6 +116,9 @@ static KEYWORD precommands[] = { ,{"sortreallocate", DoPreSortReallocate , 0, 0} #ifdef WITHFLOAT ,{"startfloat" , DoStartFloat , 0, 0} +#endif +#ifdef WITHFLINT + ,{"startpadic" , DoStartPadic , 0, 0} #endif ,{"switch" , DoPreSwitch , 0, 0} ,{"system" , DoSystem , 0, 0} @@ -7769,6 +7775,12 @@ int DoStartFloat(UBYTE *s) MesPrint("@Simultaneous use of floating point and modulus arithmetic makes no sense."); error = 1; } +#ifdef WITHFLINT + if ( AC.activePadic ) { + MesPrint("@Simultaneous use of float_ and padic_ is not allowed."); + error = 1; + } +#endif if ( AT.aux_ ) { // First, we clean up any previous floating point system. ClearfFloat(); ClearMZVTables(); @@ -7871,5 +7883,98 @@ int DoEndFloat(UBYTE *s) #endif /* #] DoEndFloat : + #[ DoStartPadic : +*/ +#ifdef WITHFLINT + +int DoStartPadic(UBYTE *s) +{ + GETIDENTITY + int error = 0; + LONG N = 0; + UBYTE *ss, *p, *pstop, c; + + if ( AP.PreSwitchModes[AP.PreSwitchLevel] != EXECUTINGPRESWITCH ) return(0); + if ( AP.PreIfStack[AP.PreIfLevel] != EXECUTINGIF ) return(0); + + if ( AR.PolyFun != 0 ) { + MesPrint("@Simultaneous use of Poly(Rat)Fun and padic_ is not allowed."); + error = 1; + } + if ( AC.ncmod != 0 ) { + MesPrint("@Simultaneous use of p-adic and modulus arithmetic makes no sense."); + error = 1; + } +#ifdef WITHFLOAT + if ( AT.aux_ != 0 ) { + MesPrint("@Simultaneous use of float_ and padic_ is not allowed."); + error = 1; + } +#endif + + while ( *s == ',' || *s == ' ' || *s == '\t' ) s++; + ss = s; + // The first parameter is the prime number + if ( FG.cTable[*s] == 1 ) { + p = s; + while ( FG.cTable[*s] == 1 ) s++; + pstop = s; + + while ( *s == ',' || *s == ' ' || *s == '\t' ) s++; + // The second parameter is the p-adic precision. + // It is signed: absolute precisions below O(p^0) is fine. + if ( *s == '+' || *s == '-' || FG.cTable[*s] == 1 ) { + UBYTE *digits = s; + while ( *digits == '+' || *digits == '-' ) digits++; + if ( FG.cTable[*digits] != 1 ) goto IllPar; + ParseSignedNumber(N,s) + } + else goto IllPar; + while ( *s == ' ' || *s == '\t' ) s++; + if ( *s != 0 ) goto IllPar; + } + else { +IllPar: + MesPrint("@Illegal parameter in %#StartPadic: %s ",ss); + error = 1; + }; + + if ( error == 0 ) { + /* Split the existing input buffer so FLINT can parse just the prime. */ + c = *pstop; + *pstop = 0; + if ( StartPadicSystem(p,N) ) error = 1; + else AO.PadicPrint = 1; + *pstop = c; + } + return(error); +} + +#endif +/* + #] DoStartPadic : + #[ DoEndPadic : +*/ +#ifdef WITHFLINT + +int DoEndPadic(UBYTE *s) +{ + int error = 0; + if ( AP.PreSwitchModes[AP.PreSwitchLevel] != EXECUTINGPRESWITCH ) return(0); + if ( AP.PreIfStack[AP.PreIfLevel] != EXECUTINGIF ) return(0); + while ( *s == ',' || *s == ' ' || *s == '\t' ) s++; + if ( *s != 0 ) { + MesPrint("@Illegal parameter in %#EndPadic instruction: %s ",s); + error = 1; + } + if ( error == 0 ) { + ClearPadicSystem(); + } + return(error); +} + +#endif +/* + #] DoEndPadic : # ] PreProcessor : */ diff --git a/sources/proces.c b/sources/proces.c index 06df6280..28fba033 100644 --- a/sources/proces.c +++ b/sources/proces.c @@ -3984,6 +3984,20 @@ SkipCount: level++; AT.WorkPointer = term + *term; if ( Chop(BHEAD term,level) ) goto GenCall; goto Return0; +#endif +#ifdef WITHFLINT + case TYPETOPADIC: + AT.WorkPointer = term + *term; + if ( ToPadic(BHEAD term,level) ) goto GenCall; + goto Return0; + case TYPEPADICTORAT: + AT.WorkPointer = term + *term; + if ( PadicToRat(BHEAD term,level) ) goto GenCall; + goto Return0; + case TYPEFROMPADIC: + AT.WorkPointer = term + *term; + if ( FromPadic(BHEAD term,level,C->lhs[level][2]) ) goto GenCall; + goto Return0; #endif } goto SkipCount; diff --git a/sources/sch.c b/sources/sch.c index 7db6efd1..1d6b8afa 100644 --- a/sources/sch.c +++ b/sources/sch.c @@ -1935,6 +1935,11 @@ int WriteInnerTerm(WORD *term, WORD first) WORD *t, *s, *s1, *s2, n, i, pow; #ifdef WITHFLOAT int FloatChars = 0; +#endif +#ifdef WITHFLINT + int PadicChars = 0; +#endif +#ifdef WITHFLOAT GETIDENTITY #endif t = term; @@ -2012,6 +2017,27 @@ int WriteInnerTerm(WORD *term, WORD first) } if ( ss >= t ) first = 1; } +#endif +#ifdef WITHFLINT +/* + Check whether there is a proper padic_ function and no raw mode. + If so, print as p-adic series. + Raw mode is indicated as AO.PadicPrint == 0. +*/ + else if ( AO.PadicPrint && AC.activePadic ) { + WORD *ss = s; + while ( ss < t ) { + if ( *ss == PADICFUN ) { + if ( ( PadicChars = PrintPadic(ss) ) != 0 ) { + TokenToLine(AO.padicspace); + first = 0; + } + break; + } + ss += ss[1]; + } + if ( ss >= t ) first = 1; + } #endif else first = 1; while ( s < t ) { @@ -2062,6 +2088,11 @@ int WriteInnerTerm(WORD *term, WORD first) if ( *s == FLOATFUN && AO.FloatPrec >= 0 && AT.aux_ != 0 ) { } else +#endif +#ifdef WITHFLINT + if ( *s == PADICFUN && AO.PadicPrint && AC.activePadic ) { + } + else #endif { if ( *s >= FUNCTION && AC.funpowers > 0 diff --git a/sources/sort.c b/sources/sort.c index 88199825..647a51a5 100644 --- a/sources/sort.c +++ b/sources/sort.c @@ -1749,6 +1749,9 @@ int AddCoef(PHEAD WORD **ps1, WORD **ps2) UWORD *OutCoef; #ifdef WITHFLOAT if ( AT.SortFloatMode ) return(AddWithFloat(BHEAD ps1,ps2)); +#endif +#ifdef WITHFLINT + if ( AT.SortPadicMode ) return(AddWithPadic(BHEAD ps1,ps2)); #endif OutCoef = AN.SoScratC; s1 = *ps1; s2 = *ps2; @@ -2362,6 +2365,9 @@ WORD Compare1(PHEAD WORD *term1, WORD *term2, WORD level) else { localPoly = 0; } #ifdef WITHFLOAT AT.SortFloatMode = 0; +#endif +#ifdef WITHFLINT + AT.SortPadicMode = 0; #endif prevorder = 0; GETSTOP(term1,s1); @@ -2692,6 +2698,18 @@ WORD Compare1(PHEAD WORD *term1, WORD *term2, WORD level) } else if ( TestFloat(s2-FUNHEAD) ) { return(-1); } } +#endif +#ifdef WITHFLINT + if ( level == 0 && c1 == PADICFUN && t1 == stopper1 && t2 == stopper2 && AC.activePadic ) { +/* + We have two PADICFUN's. Test whether they are 'legal'. +*/ + if ( TestPadic(s1-FUNHEAD) ) { + if ( TestPadic(s2-FUNHEAD) ) { AT.SortPadicMode = 3; return(0); } + else { return(1); } + } + else if ( TestPadic(s2-FUNHEAD) ) { return(-1); } + } #endif while ( s1 < t1 ) { /* @@ -2762,6 +2780,16 @@ WORD Compare1(PHEAD WORD *term1, WORD *term2, WORD level) && TestFloat(t2) && AT.aux_ != 0 ) { AT.SortFloatMode = 2; return(0); } +#endif +#ifdef WITHFLINT + if ( level == 0 && t1 < stopper1 && *t1 == PADICFUN && t1+t1[1] == stopper1 + && TestPadic(t1) && AC.activePadic ) { + AT.SortPadicMode = 1; return(0); + } + else if ( level == 0 && t2 < stopper2 && *t2 == PADICFUN && t2+t2[1] == stopper2 + && TestPadic(t2) && AC.activePadic ) { + AT.SortPadicMode = 2; return(0); + } #endif { if ( AR.SortType != SORTLOWFIRST ) { @@ -3931,6 +3959,16 @@ int MergePatches(WORD par) goto cancelled; poin[S->tree[i]] = term1; } +#endif +#ifdef WITHFLINT + else if ( AT.SortPadicMode ) { + WORD *term1, *term2; + term1 = poin[S->tree[i]]; + term2 = poin[k]; + if ( MergeWithPadic(BHEAD &term1,&term2) == 0 ) + goto cancelled; + poin[S->tree[i]] = term1; + } #endif else { r1 = *( m1 += l1 - 1 ); diff --git a/sources/startup.c b/sources/startup.c index 6f85902d..05c25117 100644 --- a/sources/startup.c +++ b/sources/startup.c @@ -1348,6 +1348,12 @@ void StartVariables(void) #ifdef WITHFLOAT AC.MaxWeight = AM.gMaxWeight = AM.ggMaxWeight = MAXWEIGHT; AC.DefaultPrecision = AM.gDefaultPrecision = AM.ggDefaultPrecision = DEFAULTPRECISION; +#endif +#ifdef WITHFLINT + AC.activePadicPrecision = 0; + AC.activePadic = 0; + AC.activePadicContext = 0; + AO.padicncoeffs = 0; #endif AC.CommuteInSet = 0; @@ -2118,6 +2124,7 @@ backtrace_fallback: ; This keeps valgrind happy. */ #ifdef WITHFLINT + ClearPadicSystem(); flint_final_cleanup_master(); #endif CleanUp(errorcode); diff --git a/sources/structs.h b/sources/structs.h index f3faeb5e..15a81dbd 100644 --- a/sources/structs.h +++ b/sources/structs.h @@ -1780,6 +1780,11 @@ struct C_const { LONG MaxWeight; /* (C) Maximum weight for MZV or Euler */ LONG tDefaultPrecision; /* (C) Default precision in bits for float_ */ LONG tMaxWeight; /* (C) Maximum weight for MZV or Euler */ +#endif +#ifdef WITHFLINT + LONG activePadicPrecision; /* (C) Precision for padic_ system */ + int activePadic; /* (C) p-adic system currently active */ + void *activePadicContext; /* (C) Active FLINT p-adic context */ #endif int cbufnum; /**< Current compiler buffer */ int AutoDeclareFlag; /** (C) Mode of looking for names. Set to NOAUTO (=0) or @@ -2110,6 +2115,9 @@ struct T_const { void *mpf_tab2; void *aux_; void *auxr_; +#endif +#ifdef WITHFLINT + void *padic_aux_; #endif NORMDATA **NormData; LONG NormDataSize; @@ -2184,6 +2192,10 @@ struct T_const { WORD FloatPos; WORD SortFloatMode; #endif +#ifdef WITHFLINT + WORD PadicPos; + WORD SortPadicMode; +#endif }; /* #] T : @@ -2379,6 +2391,11 @@ struct O_const { UBYTE *floatspace; LONG floatsize; #endif +#ifdef WITHFLINT + UBYTE *padicspace; + LONG padicsize; + LONG padicncoeffs; +#endif /*----Leave NumInBrack as first non-pointer. This is used by the checkpoints--*/ LONG NumInBrack; /* (O) For typing [] option in print */ LONG wlen; /* (O) Used to store files. */ @@ -2406,6 +2423,9 @@ struct O_const { int IndentSpace; /* For indentation in output */ #ifdef WITHFLOAT int FloatPrec; +#endif +#ifdef WITHFLINT + int PadicPrint; #endif WORD schemenum; /* for feeding a Horner scheme to Optimize */ WORD transFlag; /* () >0 indicates that translations have to be done */ diff --git a/sources/threads.c b/sources/threads.c index ebf7c8fd..29846852 100644 --- a/sources/threads.c +++ b/sources/threads.c @@ -3907,6 +3907,16 @@ int MasterMerge(void) goto cancelled; poin[S->tree[i]] = term1; } +#endif +#ifdef WITHFLINT + else if ( AT.SortPadicMode ) { + WORD *term1, *term2; + term1 = poin[S->tree[i]]; + term2 = poin[k]; + if ( MergeWithPadic(B0, &term1,&term2) == 0 ) + goto cancelled; + poin[S->tree[i]] = term1; + } #endif else { r1 = *( m1 += l1 - 1 ); @@ -4560,6 +4570,13 @@ int SortBotMerge(PHEAD0) *term1 = fun1-term1; TermFree(fun3,"MasterMerge"); } +#endif +#ifdef WITHFLINT + else if ( AT.SortPadicMode ) { + WORD *termx = term1, *termy = term2; + if ( MergeWithPadic(BHEAD &termx,&termy) == 0 ) goto cancelled; + term1 = termx; + } #endif else { r1 = *( m1 += l1 - 1 ); diff --git a/sources/transform.c b/sources/transform.c index 58a8aab4..28185f62 100644 --- a/sources/transform.c +++ b/sources/transform.c @@ -124,6 +124,12 @@ int CoTransform(UBYTE *in) MesPrint("&Illegal use of a transform statement and float_"); if ( error == 0 ) error = 1; } +#endif +#ifdef WITHFLINT + if ( (number+FUNCTION) == PADICFUN ) { + MesPrint("&Illegal use of a transform statement and padic_"); + if ( error == 0 ) error = 1; + } #endif number += MAXVARIABLES + FUNCTION; } else if ( type != CSET ) { @@ -156,6 +162,19 @@ int CoTransform(UBYTE *in) if ( error == 0 ) 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 a transform statement and padic_"); + if ( error == 0 ) error = 1; + } + } + } #endif } } @@ -765,6 +784,9 @@ int RunTransform(PHEAD WORD *term, WORD *params) hit:; #ifdef WITHFLOAT if ( *t == FLOATFUN ) goto next; +#endif +#ifdef WITHFLINT + if ( *t == PADICFUN ) goto next; #endif while ( in < t ) *out++ = *in++; tt = t + t[1]; fun = out;