From 11dd2d2d75bb6decc5fdd47699e2860745313331 Mon Sep 17 00:00:00 2001 From: cbmarini Date: Wed, 17 Jun 2026 17:48:18 +0200 Subject: [PATCH] feat: add support to work over p-adic numbers. --- check/features.frm | 1671 +++++++++++++++++++++++++++++++++++ configure.ac | 4 +- sources/Makefile.am | 3 +- sources/argument.c | 8 + sources/compcomm.c | 73 ++ sources/compiler.c | 9 + sources/declare.h | 30 +- sources/execute.c | 3 + sources/flintwrap.cc | 16 + sources/ftypes.h | 13 +- sources/function.c | 4 + sources/inivar.h | 1 + sources/normal.c | 123 +++ sources/padic.c | 2006 ++++++++++++++++++++++++++++++++++++++++++ sources/pre.c | 105 +++ sources/proces.c | 14 + sources/reshuf.c | 3 + sources/sch.c | 43 + sources/sort.c | 48 + sources/startup.c | 7 + sources/store.c | 70 +- sources/structs.h | 20 + sources/threads.c | 9 +- sources/tools.c | 5 + sources/transform.c | 22 + 25 files changed, 4291 insertions(+), 19 deletions(-) create mode 100644 sources/padic.c diff --git a/check/features.frm b/check/features.frm index 55f8317e0..ee40cb503 100644 --- a/check/features.frm +++ b/check/features.frm @@ -2382,6 +2382,1677 @@ Evaluate euler_; #pend_if wordsize == 2 assert runtime_error?("Divergent Euler sum in CalculateEuler") *--#] mzv_error_5 : +*--#[ padic_literal_units : +#StartPadic 7,20 +Off Statistics; +Format PadicPrint off; +Symbols x; +CFunctions f; +Local SYMBOLIC = padic_(7,0,20,x); +Local BARE = padic_(7,0,20,f); +Local FUNCTION = padic_(7,0,20,f(1)); +Local FRACTION = padic_(7,0,20,3/2); +Local SHORT = padic_(7,0,20,13); +Local LARGE = padic_(7,0,30,18446744073709551629); +Print; +.end +#require flint? +assert succeeded? +assert result("SYMBOLIC") =~ expr("padic_(7,0,20,x)") +assert result("BARE") =~ expr("padic_(7,0,20,f)") +assert result("FUNCTION") =~ expr("padic_(7,0,20,f(1))") +assert result("FRACTION") =~ expr("padic_(7,0,20,3/2)") +assert result("SHORT") =~ expr("padic_(7,0,20,13)") +assert result("LARGE") =~ expr("padic_(7,0,30,18446744073709551629)") +*--#] padic_literal_units : +*--#[ padic_literal_normalization : +* Signed units, prime factors and the precision cutoff. +#StartPadic 7,20 +Format PadicPrint off; +Local SMALLNEGATIVE = padic_(7,0,3,-1); +Local SMALLFACTOR = padic_(7,0,3,14); +Local SMALLZERO = padic_(7,3,3,1); +Local SMALLLARGE = padic_(7,0,12,4294967295); +Print; +.end + +* Reduction boundaries for a prime close to a power of two. +#StartPadic 2147483563,14 +Format 255; +Format PadicPrint off; +Local BOUNDBELOW = padic_(2147483563,0,14,2147483563^14-1); +Local BOUNDEQUAL = padic_(2147483563,0,14,2147483563^14); +Local BOUNDABOVE = padic_(2147483563,0,14,2147483563^14+1); +Local BOUNDFACTOR = padic_(2147483563,0,14,2147483563*(2^300+1)); +Local BOUNDTAIL = padic_(2147483563,-4,14,2147483563^18-1); +Local BOUNDLOWER = padic_(2147483563,0,14,(2^31-14*85)*2^403-1); +Local BOUNDLOWEREQUAL = padic_(2147483563,0,14,(2^31-14*85)*2^403); +Print; +.sort +Print; +.end + +* Binary units use the same normalization rules. +#StartPadic 2,100 +Format 255; +Format PadicPrint off; +Local BINARYODD = padic_(2,0,100,2^99+1); +Local BINARYEVEN = padic_(2,0,100,2*(2^90+1)); +Local BINARYZERO = padic_(2,0,100,2^100); +Print; +.sort +Print; +.end + +* Literal precision must not change subsequent rational conversion. +#StartPadic 7,20 +Format 255; +Format PadicPrint off; +Local STATELITERAL = padic_(7,0,3,-1); +.sort +Local STATERATIONAL = -1; +ToPadic; +Print; +.end +#require flint? +#pend_if wordsize == 2 +assert succeeded? +assert result("SMALLNEGATIVE") =~ expr("padic_(7,0,3,342)") +assert result("SMALLFACTOR") =~ expr("padic_(7,1,3,2)") +assert result("SMALLZERO") =~ expr("0") +assert result("SMALLLARGE") =~ expr("padic_(7,0,12,4294967295)") + +2.times do |pass| + assert result("BOUNDBELOW", pass).gsub(/\\\s*/, "") =~ expr("padic_(2147483563,0,14,#{2147483563**14-1})") + assert result("BOUNDEQUAL", pass) =~ expr("0") + assert result("BOUNDABOVE", pass) =~ expr("padic_(2147483563,0,14,1)") + assert result("BOUNDFACTOR", pass).gsub(/\\\s*/, "") =~ expr("padic_(2147483563,1,14,#{2**300+1})") + assert result("BOUNDTAIL", pass).gsub(/\\\s*/, "") =~ expr("padic_(2147483563,-4,14,#{2147483563**18-1})") + assert result("BOUNDLOWER", pass).gsub(/\\\s*/, "") =~ expr("padic_(2147483563,0,14,#{(2**31-14*85)*2**403-1})") + assert result("BOUNDLOWEREQUAL", pass).gsub(/\\\s*/, "") =~ expr("padic_(2147483563,0,14,#{(2**31-14*85)*2**403})") +end + +2.times do |pass| + assert result("BINARYODD", pass) =~ expr("padic_(2,0,100,#{2**99+1})") + assert result("BINARYEVEN", pass) =~ expr("padic_(2,1,100,#{2**90+1})") + assert result("BINARYZERO", pass) =~ expr("0") +end + +assert result("STATELITERAL") =~ expr("padic_(7,0,3,342)") +assert result("STATERATIONAL") =~ expr("padic_(7,0,20,#{7**20-1})") +*--#] padic_literal_normalization : +*--#[ padic_context_lifecycle : +* Activate and renormalize existing literals. +Symbols v; +CFunctions f; +Format 255; +Format PadicPrint off; +Local ACTNEGATIVE = padic_(7,0,12,-4294967295); +Local ACTFACTOR = padic_(7,0,12,7*4294967295); +Local ACTBOUND = padic_(7,0,12,7^12+1); +Local ACTZERO = padic_(7,12,12,18446744073709551617); +Local ACTSMALL = padic_(7,0,3,14); +Local ACTWIDE = padic_(7,-3,30,7^33-1); +.sort +#StartPadic 7,20 +Format PadicPrint off; +Print; +.sort +#EndPadic +FromPadic,f; +id f(7,v?,12,?a) = padic_(7,v,2,?a); +.sort +#StartPadic 7,20 +Format PadicPrint off; +Print; +.end + +* Replace an active context, including its multiword prime. +#StartPadic 170141183460469231731687303715884105727,2 +Format 255; +Format PadicPrint off; +Local RESTARTF = 31; +ToPadic; +Print RESTARTF; +.sort +Drop RESTARTF; +.sort +#StartPadic 7,10 +Format PadicPrint off; +Local RESTARTG = 31; +ToPadic; +Print RESTARTG; +.sort +Drop RESTARTG; +.sort +#StartPadic 170141183460469231731687303715884105727,2 +Format PadicPrint off; +Local RESTARTH = 32; +ToPadic; +Print RESTARTH; +.end +#require flint? +#pend_if wordsize == 2 +assert succeeded? +assert result("ACTNEGATIVE", 0) =~ expr("padic_(7,0,12,9546319906)") +assert result("ACTFACTOR", 0) =~ expr("padic_(7,1,12,340313809)") +assert result("ACTBOUND", 0) =~ expr("padic_(7,0,12,1)") +assert result("ACTZERO", 0) =~ expr("0") +assert result("ACTSMALL", 0) =~ expr("padic_(7,1,3,2)") +assert result("ACTWIDE", 0) =~ expr("padic_(7,-3,30,#{7**33-1})") +assert result("ACTNEGATIVE", 1) =~ expr("padic_(7,0,2,11)") +assert result("ACTFACTOR", 1) =~ expr("padic_(7,1,2,3)") +assert result("ACTBOUND", 1) =~ expr("padic_(7,0,2,1)") + +assert result("RESTARTF") =~ expr("padic_(170141183460469231731687303715884105727,0,2,31)") +assert result("RESTARTG") =~ expr("padic_(7,0,10,31)") +assert result("RESTARTH") =~ expr("padic_(170141183460469231731687303715884105727,0,2,32)") +*--#] padic_context_lifecycle : +*--#[ padic_coefficient_placement : +#- +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 + +* Active p-adic coefficients follow all monomial objects. +Format PadicPrint off; +Print; +.sort + +#EndPadic +* Inactive p-adic objects return to ordinary function ordering. +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_(7,0,10,31)") +assert result("G") =~ expr("A1*A2*cos_*cos_(pi_)*cos_(3*pi_)*mzv_*mzv_(2)*mzv_(2,1)*padic_(7,0,10,31)* + f1*f2*q*p1.p2*x1*x2") +*--#] padic_coefficient_placement : +*--#[ padic_printing : +* Series printing of signed powers and arithmetic results. +#StartPadic 7,10 +Local SERIESF = 101/7; +Local SERIESG = 7/13; +ToPadic; +.sort + +Local SERIESMUL = (10+SERIESF)*SERIESG; +Local SERIESDIV = (10+SERIESF)/SERIESG; +Print; +.end + +* Malformed factors retain ordinary printing, including Print +s. +#StartPadic 7,10 +Off Statistics; +Symbols x; +CFunctions f; +Local ORDINARYF = padic_(0,10); +Local ORDINARYNEGATIVE = -padic_(0,10); +Local ORDINARYSUM = 1+padic_(0,10); +Local ORDINARYPRODUCT = x*padic_(0,10); +Local ORDINARYCOEFFICIENT = 2*padic_(0,10); +Local ORDINARYNESTED = f(padic_(0,10)); +Local ORDINARYTYPES = padic_(0,10,x); +Local ORDINARYFRACTION = padic_(0,10,1/2); +Local ORDINARYEXTRA = padic_(0,10,1,2); +Local ORDINARYBARE = padic_; +Local ORDINARYMULTIPLE = padic_(0,10)*padic_(0,10,x); +Local ORDINARYMIXED = padic_(0,10)*padic_(7,0,10,31); +Print; +.sort +Print +s ORDINARYF,ORDINARYNEGATIVE,ORDINARYSUM,ORDINARYPRODUCT; +.end + +* Malformed prime, valuation, precision and unit arguments remain ordinary functions. +#StartPadic 7,10 +Symbols x; +Local ARGSYMBOLIC = padic_(x,0,10,1); +Local ARGFRACTION = padic_(7/2,0,10,1); +Local ARGZERO = padic_(0,0,10,1); +Local ARGONE = padic_(1,0,10,1); +Local ARGNEGATIVE = padic_(-7,0,10,1); +Local ARGLONGNEGATIVE = padic_(-170141183460469231731687303715884105727,0,10,1); +Local ARGVALUATION = padic_(7,x,10,1); +Local ARGPRECISION = padic_(7,0,x,1); +Local ARGUNIT = padic_(7,0,10,x); +Local ARGEXTRA = padic_(7,0,10,1,2); +Local ARGLEGACY = padic_(0,10,31); +Print; +.end +#require flint? +#pend_if wordsize == 2 +assert succeeded? +assert result("SERIESF") =~ expr("(3*7^-1 + 2*7^1)") +assert result("SERIESG") =~ 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("SERIESMUL") =~ 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("SERIESDIV") =~ expr("(4*7^-2 + 2*7^-1 + 3 + 6*7^1)") + +assert result("ORDINARYF", 0) =~ expr("padic_(0,10)") +assert result("ORDINARYNEGATIVE", 0) =~ expr("-padic_(0,10)") +assert result("ORDINARYSUM", 0) =~ expr("1+padic_(0,10)") +assert result("ORDINARYPRODUCT", 0) =~ expr("padic_(0,10)*x") +assert result("ORDINARYCOEFFICIENT") =~ expr("2*padic_(0,10)") +assert result("ORDINARYNESTED") =~ expr("f(padic_(0,10))") +assert result("ORDINARYTYPES") =~ expr("padic_(0,10,x)") +assert result("ORDINARYFRACTION") =~ expr("padic_(0,10,1/2)") +assert result("ORDINARYEXTRA") =~ expr("padic_(0,10,1,2)") +assert result("ORDINARYBARE") =~ expr("padic_") +assert result("ORDINARYMULTIPLE") =~ expr("padic_(0,10)*padic_(0,10,x)") +assert result("ORDINARYMIXED") =~ expr("(3+4*7^1)*padic_(0,10)") +assert result("ORDINARYF", 1) =~ expr("+padic_(0,10)") +assert result("ORDINARYNEGATIVE", 1) =~ expr("-padic_(0,10)") +assert result("ORDINARYSUM", 1) =~ expr("+1+padic_(0,10)") +assert result("ORDINARYPRODUCT", 1) =~ expr("+padic_(0,10)*x") + +assert result("ARGSYMBOLIC") =~ expr("padic_(x,0,10,1)") +assert result("ARGFRACTION") =~ expr("padic_(7/2,0,10,1)") +assert result("ARGZERO") =~ expr("padic_(0,0,10,1)") +assert result("ARGONE") =~ expr("padic_(1,0,10,1)") +assert result("ARGNEGATIVE") =~ expr("padic_(-7,0,10,1)") +assert result("ARGLONGNEGATIVE") =~ expr("padic_(-170141183460469231731687303715884105727,0,10,1)") +assert result("ARGVALUATION") =~ expr("padic_(7,x,10,1)") +assert result("ARGPRECISION") =~ expr("padic_(7,0,x,1)") +assert result("ARGUNIT") =~ expr("padic_(7,0,10,x)") +assert result("ARGEXTRA") =~ expr("padic_(7,0,10,1,2)") +assert result("ARGLEGACY") =~ expr("padic_(0,10,31)") +*--#] padic_printing : +*--#[ padic_to_rational : +* Reconstruct coefficients with positive and negative valuations. +Symbol x,j; +#StartPadic 7,5 +Local SERIESF = 101/13*sum_(j,-5,5,x^j*7^j); +ToPadic; +.sort + +PadicToRat; +Print; +.end + +* Reconstruct with negative absolute precision. +#StartPadic 7,-1 +Format PadicPrint off; +Local NEGF = 1/49; +ToPadic; +Print; +.sort +PadicToRat; +Print; +.end +#require flint? +#pend_if wordsize == 2 +assert succeeded? +assert result("SERIESF") =~ 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") + +assert result("NEGF", 0) =~ expr("padic_(7,-2,-1,1)") +assert result("NEGF", 1) =~ expr("1/49") +*--#] padic_to_rational : +*--#[ padic_laurent_reconstruction : +#- +#$prime = 175519; +#define EP "`$prime'" + +Off Statistics; +Symbol eps,x1,x2,p,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(14,4):base=`$prime'; + id padic(p?,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_laurent_reconstruction : +*--#[ padic_from_padic : +* Expand ordinary and compact units at the current execution level. +Format 255; +Format PadicPrint off; +CFunctions f,g; +* The low limb of this ordinary integer equals the PADICLIMBS type number. +Local UNITORDINARY = padic_(7,0,100,2^32+13); +Local UNITNEGATIVE = padic_(7,0,100,-(2^32+13)); +Local UNITSHORT = padic_(7,0,100,-2147483648); +Local UNITZERO = padic_(7,0,100,0); +Local UNITNONINTEGER = padic_(7,0,100,1/2); +FromPadic,f; +Print; +.sort +#StartPadic 7,100 +Local UNITCOMPACT = 2^64+13; +ToPadic; +.sort +#EndPadic +Local UNITARGUMENT = g(UNITCOMPACT); +FromPadic,f; +Print UNITCOMPACT,UNITARGUMENT; +.end + +* Conversion is independent of the active prime and context. +CFunctions f; +#StartPadic 7,10 +Local INERTF = 31; +ToPadic; +.sort +#EndPadic +Local INERTOTHER = padic_(11,0,10,1); +Local INERTLEGACY = padic_(0,10,31); +FromPadic, f; +Print; +.end +#require flint? +#pend_if wordsize == 2 +assert succeeded? +assert result("UNITORDINARY", 0) =~ expr("f(7,0,100,4294967309)") +assert result("UNITNEGATIVE", 0) =~ expr("f(7,0,100,-4294967309)") +assert result("UNITSHORT", 0) =~ expr("f(7,0,100,-2147483648)") +assert result("UNITZERO", 0) =~ expr("f(7,0,100,0)") +assert result("UNITNONINTEGER", 0) =~ expr("padic_(7,0,100,1/2)") +assert result("UNITCOMPACT") =~ expr("f(7,0,100,18446744073709551629)") +assert result("UNITARGUMENT") =~ expr("g(padic_(7,0,100,18446744073709551629))") + +assert result("INERTF") =~ expr("f(7,0,10,31)") +assert result("INERTOTHER") =~ expr("f(11,0,10,1)") +assert result("INERTLEGACY") =~ expr("padic_(0,10,31)") +*--#] padic_from_padic : +*--#[ padic_integer_conversion : +* Odd and even word counts around the integer-import boundary. +#: MaxTermSize 8000 +#StartPadic 2147483647,540 +Format 255; +Format PadicPrint off; +* Odd and even word counts below, at and above 128 FLINT limbs. +Local TRANSFERODDLOW = 2^8064+17; +Local TRANSFEREVENLOW = 2^8096+17; +Local TRANSFERODDLIMIT = 2^8128+17; +Local TRANSFEREVENLIMIT = 2^8160+17; +Local TRANSFERODDHIGH = 2^8192+17; +Local TRANSFEREVENHIGH = 2^8224+17; +Local TRANSFERLARGE = 2^16384+17; +Local TRANSFERNEGODD = padic_(2147483647,0,540,-(2^8128+17)); +Local TRANSFERNEGEVEN = padic_(2147483647,0,540,-(2^8160+17)); +Local TRANSFERNEGHIGH = padic_(2147483647,0,540,-(2^8192+17)); +ToPadic; +.sort +Multiply 3; +.sort +Local TRANSFERCOMBINED = TRANSFERODDLIMIT+TRANSFEREVENLIMIT; +Print; +.end + +* Large units survive arithmetic, cancellation and shrinking. +#StartPadic 7,100 +Format 255; +Format PadicPrint off; +CFunctions f; +Symbols x; +Local UNITF = 2^255+4294967295; +ToPadic; +.sort +Local UNITSUM = UNITF+1; +Local UNITPRODUCT = UNITF*padic_(7,0,100,2); +Local UNITQUOTIENT = UNITPRODUCT/padic_(7,0,100,2); +Local UNITZERO = UNITQUOTIENT-UNITF; +Local UNITSHRINK = UNITF-2^255; +.sort +FromPadic,f; +id f(7,0,100,x?) = x; +Print; +.end +#require flint? +#pend_if wordsize == 2 +assert succeeded? +prime = 2147483647 +precision = 540 +expected = { + "TRANSFERODDLOW" => 3*(2**8064+17), + "TRANSFEREVENLOW" => 3*(2**8096+17), + "TRANSFERODDLIMIT" => 3*(2**8128+17), + "TRANSFEREVENLIMIT" => 3*(2**8160+17), + "TRANSFERODDHIGH" => 3*(2**8192+17), + "TRANSFEREVENHIGH" => 3*(2**8224+17), + "TRANSFERLARGE" => 3*(2**16384+17), + "TRANSFERNEGODD" => prime**precision-3*(2**8128+17), + "TRANSFERNEGEVEN" => prime**precision-3*(2**8160+17), + "TRANSFERNEGHIGH" => prime**precision-3*(2**8192+17), + "TRANSFERCOMBINED" => 3*(2**8128+2**8160+34) +} +expected.each do |name, unit| + # Join continuation lines in the long decimal integers. + assert result(name).delete("\\") =~ expr("padic_(#{prime},0,#{precision},#{unit})") +end + +large = 2**255+4294967295 +assert result("UNITF") =~ expr(large.to_s) +assert result("UNITSUM") =~ expr((large+1).to_s) +assert result("UNITPRODUCT") =~ expr((2*large).to_s) +assert result("UNITQUOTIENT") =~ expr(large.to_s) +assert result("UNITZERO") =~ expr("0") +assert result("UNITSHRINK") =~ expr("4294967295") +*--#] padic_integer_conversion : +*--#[ padic_rational_conversion : +* Import large numerators, denominators and unequal padded lengths. +#StartPadic 7,100 +Format 255; +Format PadicPrint off; +* Numerator and denominator imports, including unequal padded lengths. +Local MULTINUM = (2^130+17)/3; +Local MULTIDEN = 3/(2^130+17); +Local MULTIBOTH = (2^132+17)/(2^130+3); +Local MULTINEGATIVE = -(2^132+17)/(2^130+3); +ToPadic; +.sort +Local MULTIPRODUCT = MULTIBOTH*(2^130+3); +Local MULTIMIXED = MULTIBOTH+(2^132+17)/(2^130+3); +Print; +.end + +* Denominators below and above the signed WORD boundary. +#StartPadic 2305843009213693951,1 +Format PadicPrint off; +Local WORDSMALL = 1/2147483647; +Local WORDLARGE = 1/2147483648; +Local WORDNEGATIVE = -1/2147483647; +Local WORDNUMERATOR = 2305843009213693950/2147483647; +ToPadic; +Print; +.end +#require flint? +#pend_if wordsize == 2 +assert succeeded? +modulus = 7**100 +expected = { + "MULTINUM" => [2**130+17, 3], + "MULTIDEN" => [3, 2**130+17], + "MULTIBOTH" => [2**132+17, 2**130+3], + "MULTINEGATIVE" => [-(2**132+17), 2**130+3], + "MULTIPRODUCT" => [2**132+17, 1], + "MULTIMIXED" => [2*(2**132+17), 2**130+3] +} +expected.each do |name, fraction| + # Compute the denominator inverse independently with Euclid's algorithm. + a, b = modulus, fraction[1] + x, y = 0, 1 + until b == 0 + quotient = a/b + a, b = b, a-quotient*b + x, y = y, x-quotient*y + end + assert_equal 1, a + unit = fraction[0]*x % modulus + assert result(name) =~ expr("padic_(7,0,100,#{unit})") +end + +assert result("WORDSMALL") =~ expr("padic_(2305843009213693951,0,1,2147483649)") +assert result("WORDLARGE") =~ expr("padic_(2305843009213693951,0,1,1073741824)") +assert result("WORDNEGATIVE") =~ expr("padic_(2305843009213693951,0,1,2305843007066210302)") +assert result("WORDNUMERATOR") =~ expr("padic_(2305843009213693951,0,1,2305843007066210302)") +*--#] padic_rational_conversion : +*--#[ padic_large_primes : +* Conversion, arithmetic and reconstruction with a multiword prime. +#StartPadic 170141183460469231731687303715884105727,2 +Format 255; +Format PadicPrint off; +Local MULTIF = 2^80; +ToPadic; +Print; +.sort +Local MULTISUM = MULTIF+1; +Local MULTIPRODUCT = MULTIF*MULTIF; +Local MULTIQUOTIENT = MULTIPRODUCT/MULTIF; +Local MULTIZERO = MULTIF-MULTIF; +Print MULTISUM,MULTIPRODUCT,MULTIQUOTIENT,MULTIZERO; +.sort +Skip; +NSkip MULTIF; +PadicToRat; +Print MULTIF; +.end + +* A prime immediately above the signed WORD range. +#StartPadic 2147483659,2 +Format PadicPrint off; +Local WORDF = 31; +ToPadic; +Print; +.sort +Local WORDG = WORDF+1; +Print WORDG; +.end +#require flint? +#pend_if wordsize == 2 +assert succeeded? +assert result("MULTIF", 0) =~ expr("padic_(170141183460469231731687303715884105727,0,2,1208925819614629174706176)") +assert result("MULTISUM") =~ expr("padic_(170141183460469231731687303715884105727,0,2,1208925819614629174706177)") +assert result("MULTIPRODUCT") =~ expr("padic_(170141183460469231731687303715884105727,0,2,1461501637330902918203684832716283019655932542976)") +assert result("MULTIQUOTIENT") =~ expr("padic_(170141183460469231731687303715884105727,0,2,1208925819614629174706176)") +assert result("MULTIZERO") =~ expr("0") +assert result("MULTIF", 1) =~ expr("1208925819614629174706176") + +assert result("WORDF") =~ expr("padic_(2147483659,0,2,31)") +assert result("WORDG") =~ expr("padic_(2147483659,0,2,32)") +*--#] padic_large_primes : +*--#[ padic_rational_scaling : +#StartPadic 7,3 +Format PadicPrint off; +* Rational scaling is exact; the input's stored precision remains a bound. +Local FRACTIONUNIT = padic_(7,0,3,3); +Local FRACTIONPOS = padic_(7,1,3,3); +Local FRACTIONNEG = padic_(7,-2,3,3); +Local FRACTIONLOW = padic_(7,1,2,3); +Local FRACTIONHIGH = padic_(7,1,12,3); +Local FRACTIONZERO = padic_(7,0,3,0); +Local FRACTIONABOVE = padic_(7,4,12,3); +* Divide the represented unit exactly, including beyond the cached powers. +Local DEEPEDGE = padic_(7,-3,3,117648); +Local DEEPOUTSIDE = padic_(7,-4,3,823542); +Local DEEPLOW = padic_(7,-3,1,2400); +Local DEEPEXACT = padic_(7,-3,3,2); +* Negation preserves stored precision and complements modulo p^(N-v). +Local NEGATEUNIT = padic_(7,0,3,3); +Local NEGATEPOS = padic_(7,1,3,3); +Local NEGATENEG = padic_(7,-2,3,3); +Local NEGATEEDGE = padic_(7,-3,3,117648); +Local NEGATELOW = padic_(7,-3,1,2400); +Local NEGATEHIGH = padic_(7,-2,4,100000); +Local NEGATEOUTSIDE = padic_(7,-4,3,823542); +Local NEGATESIGNED = padic_(7,-2,3,-3); +Local NEGATEDIVISIBLE = padic_(7,-2,3,21); +* Cap wide inputs and shift negative valuations. +Local CAPHIGH = padic_(7,0,4,2000); +Local CAPHIGHNEG = padic_(7,-1,4,10000); +Local CAPUNIT = padic_(7,-2,3,3); +.sort +Local FRACTIONU = FRACTIONUNIT/2; +Local FRACTIONUM = -FRACTIONUNIT/2; +Local FRACTIONP = FRACTIONPOS/2; +Local FRACTIONN = FRACTIONNEG/2; +Local FRACTIONNM = -FRACTIONNEG/2; +Local FRACTIONL = FRACTIONLOW/2; +Local FRACTIONH = FRACTIONHIGH/2; +Local FRACTIONZ = FRACTIONZERO/2; +Local FRACTIONA = FRACTIONABOVE/2; +Local FRACTIONNUM = 7*FRACTIONUNIT/2; +Local FRACTIONDEN = FRACTIONUNIT/7; +Local FRACTIONOUTSIDE = FRACTIONUNIT/8; +Local DEEPE = DEEPEDGE/2; +Local DEEPM = -DEEPEDGE/2; +Local DEEPO = DEEPOUTSIDE/2; +Local DEEPL = DEEPLOW/2; +Local DEEPX = DEEPEXACT/2; +Local NEGATEU = -NEGATEUNIT; +Local NEGATEP = -NEGATEPOS; +Local NEGATEN = -NEGATENEG; +Local NEGATEE = -NEGATEEDGE; +Local NEGATEL = -NEGATELOW; +Local NEGATEH = -NEGATEHIGH; +Local NEGATEO = -NEGATEOUTSIDE; +Local NEGATES = -NEGATESIGNED; +Local NEGATED = -NEGATEDIVISIBLE; +Local CAPH = CAPHIGH/2; +Local CAPN = CAPHIGHNEG/2; +Local CAPZ = 343*CAPUNIT; +Print; +.end + +* Signed exact integers shift valuations and precision. +#StartPadic 7,2 +Format PadicPrint off; +* Signed integer scaling uses the full unit modulus and retains input bounds. +Local INTEGERUNIT = padic_(7,-1,2,1); +Local INTEGERPOS = padic_(7,1,2,1); +Local INTEGEREDGE = padic_(7,-2,2,48); +Local INTEGERLOW = padic_(7,-1,1,1); +Local INTEGERHIGH = padic_(7,-1,3,342); +Local INTEGEROUTSIDE = padic_(7,-3,2,1); +.sort +Local INTEGERP = 2*INTEGERUNIT; +Local INTEGERN = -2*INTEGERUNIT; +Local INTEGERF = -8*INTEGERUNIT; +Local INTEGERV = -2*INTEGERPOS; +Local INTEGERE = -2*INTEGEREDGE; +Local INTEGERL = -2*INTEGERLOW; +Local INTEGERH = -2*INTEGERHIGH; +Local INTEGERO = -2*INTEGEROUTSIDE; +Local INTEGERD = 7*INTEGERUNIT; +Local INTEGERI = INTEGERLOW; +Local INTEGERZ = 0*INTEGERUNIT; +Print INTEGERP,INTEGERN,INTEGERF,INTEGERV,INTEGERE,INTEGERL,INTEGERH,INTEGERO,INTEGERD,INTEGERI,INTEGERZ; +.end + +* Exact rational factors retain guard digits and shift absolute precision. +#StartPadic 101,12 +Format 255; +Format PadicPrint off; +* Exact rational multipliers must provide guard digits for negative valuations. +Local NEGATE = -padic_(101,-2,12,1); +Local SMALL = -1/2*padic_(101,-2,12,1); +* The required arithmetic modulus lies beyond the context's cached powers. +Local DEEP = -1/2*padic_(101,-13,12,1); +* Numerator and denominator exceed the small-rational shortcut's bounds. +Local GENERAL = -102/103*padic_(101,-2,12,1); +Local WIDE = -1/2*padic_(101,-2,20,101^21+7); +Local POSITIVE = -2/3*padic_(101,1,20,101^18+1); +Local SHORT = -1/2*padic_(101,1,3,1); +Local NUMPOWER = -101/103*padic_(101,-2,3,1); +Local LOWER = padic_(101,1,3,1)/101; +Local HIGHER = 101*padic_(101,-1,3,1); +Local NEGATIVE = padic_(101,0,3,1)/101^5; +* Test product valuations below, at, and above the active precision. +Local BELOWPRECISION = 101*padic_(101,10,20,1); +Local ATPRECISION = 101^2*padic_(101,10,20,1); +Local ABOVEPRECISION = 101^3/2*padic_(101,10,20,1); +Local NEGATIVEZERO = 101^20*padic_(101,-2,12,1); +Print; +.end +#require flint? +assert succeeded? +expected = { + "FRACTIONU" => "padic_(7,0,3,173)", + "FRACTIONUM" => "padic_(7,0,3,170)", + "FRACTIONP" => "padic_(7,1,3,26)", + "FRACTIONN" => "padic_(7,-2,3,8405)", + "FRACTIONNM" => "padic_(7,-2,3,8402)", + "FRACTIONL" => "padic_(7,1,2,5)", + "FRACTIONH" => "padic_(7,1,3,26)", + "FRACTIONZ" => "0", + "FRACTIONA" => "0", + "FRACTIONNUM" => "padic_(7,1,3,26)", + "FRACTIONDEN" => "padic_(7,-1,2,3)", + "FRACTIONOUTSIDE" => "padic_(7,0,3,129)", +} +expected.each do |name, value| + assert result(name) =~ expr(value) +end + +assert result("DEEPE") =~ expr("padic_(7,-3,3,58824)") +assert result("DEEPM") =~ expr("padic_(7,-3,3,58825)") +assert result("DEEPO") =~ expr("padic_(7,-4,3,411771)") +assert result("DEEPL") =~ expr("padic_(7,-3,1,1200)") +assert result("DEEPX") =~ expr("padic_(7,-3,3,1)") + +expected = { + "NEGATEU" => "padic_(7,0,3,340)", + "NEGATEP" => "padic_(7,1,3,46)", + "NEGATEN" => "padic_(7,-2,3,16804)", + "NEGATEE" => "padic_(7,-3,3,1)", + "NEGATEL" => "padic_(7,-3,1,1)", + "NEGATEH" => "padic_(7,-2,3,842)", + "NEGATEO" => "padic_(7,-4,3,1)", + "NEGATES" => "padic_(7,-2,3,3)", + "NEGATED" => "padic_(7,-1,3,2398)" +} +expected.each do |name, value| + assert result(name) =~ expr(value) +end + +expected = { + "INTEGERP" => "padic_(7,-1,2,2)", + "INTEGERN" => "padic_(7,-1,2,341)", + "INTEGERF" => "padic_(7,-1,2,335)", + "INTEGERV" => "padic_(7,1,2,5)", + "INTEGERE" => "padic_(7,-2,2,2305)", + "INTEGERL" => "padic_(7,-1,1,47)", + "INTEGERH" => "padic_(7,-1,2,2)", + "INTEGERO" => "padic_(7,-3,2,16805)", + "INTEGERD" => "padic_(7,0,2,1)", + "INTEGERI" => "padic_(7,-1,1,1)", + "INTEGERZ" => "0" +} +expected.each do |name, value| + assert result(name) =~ expr(value) +end + +expected = { + "CAPH" => "padic_(7,0,3,314)", + "CAPN" => "padic_(7,-1,3,198)", + "CAPZ" => "padic_(7,1,3,3)", +} +expected.each do |name, value| + assert result(name) =~ expr(value) +end + +expected = { + "NEGATE" => [-2, 12, -1, 1], + "SMALL" => [-2, 12, -1, 2], + "DEEP" => [-13, 12, -1, 2], + "GENERAL" => [-2, 12, -102, 103], + "WIDE" => [-2, 12, -7, 2], + "POSITIVE" => [1, 12, -2, 3], + "SHORT" => [1, 3, -1, 2], + "NUMPOWER" => [-1, 4, -1, 103], + "LOWER" => [0, 2, 1, 1], + "HIGHER" => [0, 4, 1, 1], + "NEGATIVE" => [-5, -2, 1, 1], + "BELOWPRECISION" => [11, 12, 1, 1], +} +expected.each do |name, value| + valuation, precision, numerator, denominator = value + modulus = 101**(precision-valuation) + # Compute the exact multiplier's residue independently with Euclid. + a, b = modulus, denominator + x, y = 0, 1 + until b == 0 + quotient = a/b + a, b = b, a-quotient*b + x, y = y, x-quotient*y + end + assert_equal 1, a + unit = numerator*x % modulus + assert result(name).gsub(/\\\s*/, "") =~ expr("padic_(101,#{valuation},#{precision},#{unit})") +end +assert result("ATPRECISION") =~ expr("0") +assert result("ABOVEPRECISION") =~ expr("0") +assert result("NEGATIVEZERO") =~ expr("0") +*--#] padic_rational_scaling : +*--#[ padic_addition : +#StartPadic 7,3 +Format PadicPrint off; +* Unequal valuations do not restore digits absent from the lower-precision input. +Local ALIGNA = padic_(7,-2,3,1)+padic_(7,-1,3,2); +Local ALIGNB = padic_(7,-1,3,2)+padic_(7,-2,3,1); +Local ALIGNC = padic_(7,0,3,342)+padic_(7,1,3,48); +Local ALIGNE = padic_(7,-2,3,16806)+padic_(7,-1,3,2400); +Local ALIGNSHIFT = padic_(7,-2,3,1)+padic_(7,1,3,1); +Local ALIGNMISS = padic_(7,-2,3,1)+padic_(7,2,3,1); +Local ALIGNOUTSIDE = padic_(7,-4,3,1)+padic_(7,-3,3,1); +Local ALIGNLOW = padic_(7,-2,1,1)+padic_(7,-1,3,2); +Local ALIGNHIGH = padic_(7,-2,4,100000)+padic_(7,-1,3,2); +Local ALIGNZ = padic_(7,0,3,0)+padic_(7,-1,3,2); +Local ALIGNCANCEL = padic_(7,-1,3,1)+padic_(7,-1,3,2400); +Local ALIGNRAISE = padic_(7,-1,3,1)+padic_(7,-1,3,6); +* All sign combinations and mixed exact rational coefficients. +Local SIGNPOSITIVE = padic_(7,0,3,5)+padic_(7,0,3,3); +Local SIGNLEFT = -padic_(7,0,3,5)+padic_(7,0,3,3); +Local SIGNRIGHT = padic_(7,0,3,5)-padic_(7,0,3,3); +Local SIGNBOTH = -padic_(7,0,3,5)-padic_(7,0,3,3); +Local SIGNCANCEL = padic_(7,0,3,5)-padic_(7,0,3,5); +Local SIGNMIXED = -padic_(7,0,3,5)+1/2; +Local SIGNREVERSE = 1/2-padic_(7,0,3,5); +Local SIGNNEGATIVE = -1/2-padic_(7,0,3,5); +* Carry, wrap, cancellation and unequal stored precision. +Local REDUCECANCEL = padic_(7,-1,3,1200)+padic_(7,-1,3,1201); +Local REDUCEWRAP = padic_(7,-1,3,1200)+padic_(7,-1,3,1202); +Local REDUCEDIVISIBLE = padic_(7,-1,3,12)+padic_(7,-1,3,37); +Local REDUCEMIXED = padic_(7,-1,2,10)+padic_(7,-1,3,10); +Print; +.end + +* Addition uses the smaller input precision, including N <= 0. +#StartPadic 101,12 +Format PadicPrint off; +Local PRECADDUNIT = padic_(101,0,3,1)+padic_(101,0,5,2); +Local PRECADDVAL = padic_(101,-2,3,1)+padic_(101,1,5,2); +Local PRECADDRAT = padic_(101,0,3,1)+1; +Local PRECREVERSE = 1+padic_(101,0,3,1); +* Absolute precision may legitimately become zero or negative. +Local PRECADDZERO = padic_(101,-3,0,1)+padic_(101,-2,1,2); +Local PRECADDNEG = padic_(101,-3,-1,1)+padic_(101,-2,0,2); +Print; +.end + +* Cancellation of a converted negative-valuation coefficient. +Off Statistics; +#StartPadic 101,12 +Local SELFA = 1/101^2; +ToPadic; +.sort +Local SELFD = SELFA-SELFA; +Print SELFD; +.sort +#$n = termsin_(SELFD); +#if (`$n' != 0) + #terminate 1 +#endif +.end +#require flint? +assert succeeded? +expected = { + "ALIGNA" => "padic_(7,-2,3,15)", + "ALIGNB" => "padic_(7,-2,3,15)", + "ALIGNC" => "padic_(7,0,3,335)", + "ALIGNE" => "padic_(7,-2,3,16799)", + "ALIGNSHIFT" => "padic_(7,-2,3,344)", + "ALIGNMISS" => "padic_(7,-2,3,2402)", + "ALIGNOUTSIDE" => "padic_(7,-4,3,8)", + "ALIGNLOW" => "padic_(7,-2,1,15)", + "ALIGNHIGH" => "padic_(7,-2,3,15979)", + "ALIGNZ" => "padic_(7,-1,3,2)", + "ALIGNCANCEL" => "0", + "ALIGNRAISE" => "padic_(7,0,3,1)" +} +expected.each do |name, value| + assert result(name) =~ expr(value) +end + +expected = { + "SIGNPOSITIVE" => "padic_(7,0,3,8)", + "SIGNLEFT" => "padic_(7,0,3,341)", + "SIGNRIGHT" => "padic_(7,0,3,2)", + "SIGNBOTH" => "padic_(7,0,3,335)", + "SIGNCANCEL" => "0", + "SIGNMIXED" => "padic_(7,0,3,167)", + "SIGNREVERSE" => "padic_(7,0,3,167)", + "SIGNNEGATIVE" => "padic_(7,0,3,166)", +} +expected.each do |name, value| + assert result(name) =~ expr(value) +end + +expected = { + "REDUCECANCEL" => "0", + "REDUCEWRAP" => "padic_(7,-1,3,1)", + "REDUCEDIVISIBLE" => "padic_(7,1,3,1)", + "REDUCEMIXED" => "padic_(7,-1,2,20)" +} +expected.each do |name, value| + assert result(name) =~ expr(value) +end + +expected = { + "PRECADDUNIT" => [0, 3, 3], + "PRECADDVAL" => [-2, 3, 2060603], + "PRECADDRAT" => [0, 3, 2], + "PRECREVERSE" => [0, 3, 2], + "PRECADDZERO" => [-3, 0, 203], + "PRECADDNEG" => [-3, -1, 203], +} +expected.each do |name, value| + assert result(name) =~ expr("padic_(101,#{value.join(",")})") +end + +assert result("SELFD") =~ expr("0") +*--#] padic_addition : +*--#[ padic_multiplication : +* Product precision is bounded by each input error shifted by the other valuation. +#StartPadic 101,12 +Format PadicPrint off; +Local MULNEG = padic_(101,-2,12,1)*padic_(101,0,12,2); +Local MULLOW = padic_(101,2,5,1)*padic_(101,-1,4,2); +* Absolute precision may legitimately become zero or negative. +Local MULZERO = padic_(101,-2,2,1)*padic_(101,-1,2,2); +Local MULNEGATIVE = padic_(101,-3,1,1)*padic_(101,-2,3,2); +Print; +.end +#require flint? +assert succeeded? +expected = { + "MULNEG" => [-2, 10, 2], + "MULLOW" => [1, 4, 2], + "MULZERO" => [-3, 0, 2], + "MULNEGATIVE" => [-5, -1, 2], +} +expected.each do |name, value| + assert result(name) =~ expr("padic_(101,#{value.join(",")})") +end +*--#] padic_multiplication : +*--#[ padic_division : +* Inversion and division retain both independent precision bounds, including N <= 0. +#StartPadic 101,12 +Format PadicPrint off; +Local INV = 1/padic_(101,2,5,1); +Local INVNEG = 1/padic_(101,-2,5,1); +Local DIV = padic_(101,2,5,1)/padic_(101,1,4,1); +Local INVUNIT = 1/padic_(101,2,7,3); +Local INVNEGUNIT = 1/padic_(101,-2,5,13); +* Exercise each independent division precision bound with nontrivial units. +Local DIVNUM = padic_(101,-1,5,17)/padic_(101,2,10,13); +Local DIVDEN = padic_(101,-1,10,17)/padic_(101,2,7,13); +* Absolute precision may legitimately become zero or negative. +Local INVZERO = 1/padic_(101,2,4,1); +Local INVNEGATIVE = 1/padic_(101,3,4,1); +Local DIVZERO = padic_(101,0,2,1)/padic_(101,2,4,1); +Local DIVNEGATIVE = padic_(101,0,1,1)/padic_(101,3,5,1); +Print; +.end +#require flint? +assert succeeded? +expected = { + "INV" => [-2, 1, 1], + "INVNEG" => [2, 9, 1], + "DIV" => [1, 4, 1], + "INVZERO" => [-2, 0, 1], + "INVNEGATIVE" => [-3, -2, 1], + "DIVZERO" => [-2, 0, 1], + "DIVNEGATIVE" => [-3, -2, 1] +} +expected.each do |name, value| + assert result(name) =~ expr("padic_(101,#{value.join(",")})") +end +quotients = { + "INVUNIT" => [-2, 3, 1, 3], + "INVNEGUNIT" => [2, 9, 1, 13], + "DIVNUM" => [-3, 3, 17, 13], + "DIVDEN" => [-3, 2, 17, 13] +} +quotients.each do |name, value| + valuation, precision, numerator, denominator = value + modulus = 101**(precision-valuation) + a, b = modulus, denominator + x, y = 0, 1 + until b == 0 + quotient = a/b + a, b = b, a-quotient*b + x, y = y, x-quotient*y + end + assert_equal 1, a + unit = numerator*x % modulus + assert result(name) =~ expr("padic_(101,#{valuation},#{precision},#{unit})") +end +*--#] padic_division : +*--#[ padic_nonpositive_precision : +#StartPadic 101,0 +Format PadicPrint off; +Local SMALL = -1/2*padic_(101,-3,2,3); +Local GENERAL = -102/103*padic_(101,-3,2,3); +Local CONVERT = -1/(2*101^3); +ToPadic; +Print; +.end +#StartPadic 101,-2 +Format PadicPrint off; +Local SMALL = -1/2*padic_(101,-5,-1,3); +Local GENERAL = -102/103*padic_(101,-5,-1,3); +Local CONVERT = -1/(2*101^5); +ToPadic; +Print; +.end +#require flint? +assert succeeded? +# Contexts with N <= 0 cannot use the positive-precision shortcut. +[[-3, 0], [-5, -2]].each_with_index do |bounds, pass| + valuation, precision = bounds + modulus = 101**(precision-valuation) + { "SMALL" => [-3, 2], "GENERAL" => [-306, 103], "CONVERT" => [-1, 2] }.each do |name, fraction| + a, b = modulus, fraction[1] + x, y = 0, 1 + until b == 0 + quotient = a/b + a, b = b, a-quotient*b + x, y = y, x-quotient*y + end + assert_equal 1, a + unit = fraction[0]*x % modulus + assert result(name, pass) =~ expr("padic_(101,#{valuation},#{precision},#{unit})") + end +end +*--#] padic_nonpositive_precision : +*--#[ padic_precision_isolation : +#StartPadic 101,12 +Format 255; +Format PadicPrint off; +* Each module starts with only its new expression, reusing arithmetic scratch. +Local LOWADD = padic_(101,0,3,1)+padic_(101,0,5,2); +Print; +.sort +Drop LOWADD; +Local AFTERADD = -1/2*padic_(101,-2,12,3); +Print; +.sort +Drop AFTERADD; +Local LOWMUL = padic_(101,2,5,3)*padic_(101,-1,4,2); +Print; +.sort +Drop LOWMUL; +Local AFTERMUL = padic_(101,0,12,3)+padic_(101,0,12,4); +Print; +.sort +Drop AFTERMUL; +Local LOWDIV = padic_(101,-1,5,17)/padic_(101,2,10,13); +Print; +.sort +Drop LOWDIV; +Local AFTERDIV = padic_(101,0,12,5)*padic_(101,0,12,6); +Print; +.sort +Drop AFTERDIV; +Local LOWINV = 1/padic_(101,2,5,3); +Print; +.sort +Drop LOWINV; +Local AFTERINV = padic_(101,-1,12,5)/padic_(101,0,12,3); +Print; +.sort +Drop AFTERINV; +Local LOWRAT = padic_(101,0,3,1)/101^5; +Print; +.sort +Drop LOWRAT; +Local AFTERRAT = 1/padic_(101,0,12,3); +Print; +.sort +Drop AFTERRAT; +Local LOWCONVERT = padic_(101,0,3,1)/101^5; +.sort +Drop LOWCONVERT; +Local RESET = -1; +ToPadic; +Print RESET; +.end +#require flint? +assert succeeded? +expected = { + "LOWADD" => [0, 3, 3, 1], + "AFTERADD" => [-2, 12, -3, 2], + "LOWMUL" => [1, 4, 6, 1], + "AFTERMUL" => [0, 12, 7, 1], + "LOWDIV" => [-3, 3, 17, 13], + "AFTERDIV" => [0, 12, 30, 1], + "LOWINV" => [-2, 1, 1, 3], + "AFTERINV" => [-1, 11, 5, 3], + "LOWRAT" => [-5, -2, 1, 1], + "AFTERRAT" => [0, 12, 1, 3], + "RESET" => [0, 12, -1, 1] +} +expected.each do |name, value| + valuation, precision, numerator, denominator = value + modulus = 101**(precision-valuation) + a, b = modulus, denominator + x, y = 0, 1 + until b == 0 + quotient = a/b + a, b = b, a-quotient*b + x, y = y, x-quotient*y + end + assert_equal 1, a + unit = numerator*x % modulus + assert result(name) =~ expr("padic_(101,#{valuation},#{precision},#{unit})") +end +*--#] padic_precision_isolation : +*--#[ padic_sort_in_memory : +* Keep the default sort buffers so this identity exercises larger in-memory merges. +#- +Off Statistics; +#StartPadic 7,10 +Symbols x1,...,x4; +Local F = (x1+padic_(7,0,10,1)*x2+x3+padic_(7,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") +*--#] padic_sort_in_memory : +*--#[ padic_sort_coefficients : +* Malformed functions remain monomial factors, not coefficients. +#StartPadic 7,100 +Off Statistics; +Format 255; +Format PadicPrint off; +Symbols x; +Local FACTORBOTH = padic_(0,10)*(padic_(7,0,100,1)+padic_(7,0,100,2)); +Local FACTORRIGHT = padic_(0,10)*(1+padic_(7,0,100,2)); +Local FACTORLEFT = padic_(0,10)*(padic_(7,0,100,2)+1); +Local FACTORSYMBOLIC = padic_(7,0,10,x)*(padic_(7,0,100,1)+padic_(7,0,100,2)); +Local FACTORBARE = padic_*(padic_(7,0,100,1)+padic_(7,0,100,2)); +Local FACTORLARGE = padic_(0,10)*(padic_(7,0,100,2^130+13)+padic_(7,0,100,2^130+17)); +Print; +.end + +* Mixed coefficients have consistent ordering in both directions. +#StartPadic 7,5 +Off Statistics; +On LowFirst; +Symbols a; +Local ORDERF = 1/2*a+1; +ToPadic; +.sort +Local ORDERFORWARD = ORDERF-1/2*a-1; +Local ORDERREVERSE = -1-1/2*a+ORDERF; +Local ORDERNEGATED = 1+1/2*a-ORDERF; +Print ORDERFORWARD,ORDERREVERSE,ORDERNEGATED; +.end + +* Functions, dot products, powers and symbols share the same ordering. +#StartPadic 7,100 +Off Statistics; +On LowFirst; +Symbols a,b; +Vectors p,q; +CFunctions f; +Local OBJECTF = (2^130+13)*a+1; +Local OBJECTFUN = 1/2*f(a)+1; +Local OBJECTDOT = 1/3*p.q+1; +Local OBJECTPOWER = 1/5/a+1; +Local OBJECTMONOMIAL = 1/2*a*b+1; +ToPadic; +.sort +Local OBJECTLARGE = OBJECTF-(2^130+13)*a-1; +Local OBJECTFUNCTION = OBJECTFUN-1/2*f(a)-1; +Local OBJECTDOTPRODUCT = OBJECTDOT-1/3*p.q-1; +Local OBJECTNEGPOWER = OBJECTPOWER-1/5/a-1; +Local OBJECTPRODUCT = OBJECTMONOMIAL-1/2*a*b-1; +Print OBJECTLARGE,OBJECTFUNCTION,OBJECTDOTPRODUCT,OBJECTNEGPOWER,OBJECTPRODUCT; +.sort +On HighFirst; +Local OBJECTHIGH = OBJECTF-(2^130+13)*a-1; +Print OBJECTHIGH; +.end +#require flint? +assert succeeded? +assert result("FACTORBOTH") =~ expr("padic_(0,10)*padic_(7,0,100,3)") +assert result("FACTORRIGHT") =~ expr("padic_(0,10)*padic_(7,0,100,3)") +assert result("FACTORLEFT") =~ expr("padic_(0,10)*padic_(7,0,100,3)") +assert result("FACTORSYMBOLIC") =~ expr("padic_(7,0,10,x)*padic_(7,0,100,3)") +assert result("FACTORBARE") =~ expr("padic_*padic_(7,0,100,3)") +assert result("FACTORLARGE") =~ expr("padic_(0,10)*padic_(7,0,100,#{2**131+30})") + +assert result("ORDERFORWARD") =~ expr("0") +assert result("ORDERREVERSE") =~ expr("0") +assert result("ORDERNEGATED") =~ expr("0") + +assert result("OBJECTLARGE") =~ expr("0") +assert result("OBJECTFUNCTION") =~ expr("0") +assert result("OBJECTDOTPRODUCT") =~ expr("0") +assert result("OBJECTNEGPOWER") =~ expr("0") +assert result("OBJECTPRODUCT") =~ expr("0") +assert result("OBJECTHIGH") =~ expr("0") +*--#] padic_sort_coefficients : +*--#[ padic_sort_patch_merge : +* Force small patches so these cases exercise merging as well as in-memory sorting. +* Malformed factors with small and multiword units. +#: termsinsmall 16 +On fewerstats 1; +#StartPadic 7,100 +Format 255; +Format PadicPrint off; +Symbols x1,...,x16; +Local FACTORF = padic_(0,10)*(x1+...+x16)^4*padic_(7,0,100,3); +Local FACTORMIX = padic_(0,10)*(x1+x2)^8*(1+padic_(7,0,100,2)); +Local FACTORLARGE = padic_(0,10)*(x1+...+x16)^4*padic_(7,0,100,2^130+13); +#do i=1,16 + id x`i' = 1; +#enddo +Print; +.end + +* Mixed rational and p-adic monomials merge in LowFirst order. +#: termsinsmall 16 +#StartPadic 7,100 +Off Statistics; +On LowFirst; +Symbols a,b1,...,b8; +Local ORDERF = (1/2*a+1)*(1+b1+...+b8)^3; +ToPadic; +.sort +Local ORDERZERO = ORDERF-(1/2*a+1)*(1+b1+...+b8)^3; +Print ORDERZERO; +.end + +#require flint? +assert succeeded? +assert result("FACTORF") =~ expr("padic_(0,10)*padic_(7,0,100,196608)") +assert result("FACTORMIX") =~ expr("padic_(0,10)*padic_(7,0,100,768)") +assert result("FACTORLARGE") =~ expr("padic_(0,10)*padic_(7,0,100,#{16**4*(2**130+13)})") + +assert result("ORDERZERO") =~ expr("0") + +*--#] padic_sort_patch_merge : +*--#[ padic_sort_patch_signed_coefficients : +* Signed coefficients and exact rationals cancel across patches. +* ParFORM's master merge currently adds rational carriers instead of p-adic values. +#: termsinsmall 16 +#StartPadic 7,3 +Format PadicPrint off; +Symbols x,b1,...,b8; +Local SIGNF = (padic_(7,0,3,5)*x+padic_(7,0,3,3))*(1+b1+...+b8)^3; +.sort +Local SIGNMIXED = SIGNF-(5*x+3)*(1+b1+...+b8)^3; +Local SIGNCANCEL = SIGNF+(-padic_(7,0,3,5)*x-padic_(7,0,3,3))*(1+b1+...+b8)^3; +Local SIGNNEGATIVE = -SIGNF; +id x = 1; +#do i=1,8 + id b`i' = 1; +#enddo +Print; +.end +#require flint? +#pend_if mpi? +assert succeeded? +assert result("SIGNF") =~ expr("padic_(7,0,3,1)") +assert result("SIGNMIXED") =~ expr("0") +assert result("SIGNCANCEL") =~ expr("0") +assert result("SIGNNEGATIVE") =~ expr("padic_(7,0,3,342)") +*--#] padic_sort_patch_signed_coefficients : +*--#[ padic_sort_patch_precision : +* Patch sums retain the smaller stored precision. +* ParFORM's master merge currently adds rational carriers instead of p-adic values. +#: termsinsmall 16 +#StartPadic 101,12 +Format PadicPrint off; +Symbols x,b1,...,b8; +Local PRECBOTH = (padic_(101,0,3,1)*x+padic_(101,0,5,2))*(1+b1+...+b8)^3; +Local PRECMIXED = (padic_(101,0,3,1)*x+1)*(1+b1+...+b8)^3; +Local PRECREVERSE = (1+padic_(101,0,3,1)*x)*(1+b1+...+b8)^3; +.sort +id x = 1; +#do i=1,8 + id b`i' = 1; +#enddo +Print; +.end +#require flint? +#pend_if mpi? +assert succeeded? +assert result("PRECBOTH") =~ expr("padic_(101,0,3,2187)") +assert result("PRECMIXED") =~ expr("padic_(101,0,3,1458)") +assert result("PRECREVERSE") =~ expr("padic_(101,0,3,1458)") +*--#] padic_sort_patch_precision : +*--#[ padic_sort_patch_growth : +#: termsinsmall 16 +On fewerstats 1; +#StartPadic 7,10 +Format PadicPrint off; +Symbol x1,...,x16; + +Local F = (x1+padic_(7,0,10,1)*x2+x3+padic_(7,0,10,1)*x4)^5; +id x1 = 1-x2-x3-x4; +.sort + +Local CORNER = x1+...+x15+padic_(7,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_(7,0,10,188220513) + x15 + x14 + x13 + x12 + x11 + x10 + x9 + x8 + + x7 + x6 + x5 + x4 + x3 + x2 + x1") +*--#] padic_sort_patch_growth : +*--#[ padic_inactive_algebra : +* Storage-size coverage lives in padic_storage_size; this fold checks ordinary algebra. +* Existing coefficients and new literals remain inert after EndPadic. +CFunctions f; +#StartPadic 7,12 +Format PadicPrint off; +Global STOREDCOMPACT = padic_(7,0,12,4294967295); +.sort +#EndPadic +Global STOREDPOSITIVE = padic_(7,0,12,4294967295); +Global STOREDNEGATIVE = padic_(7,0,12,-4294967295); +Global STOREDNESTED = f(padic_(7,0,12,18446744073709551617)); +Global STOREDSMALL = padic_(7,0,12,-2147483648); +Print; +.end + +* Inactive units follow ordinary function multiplication and ordering. +CFunctions f,g; +Global ORDEROLD = padic_(7,0,12,4294967295); +.sort +Local ORDERF = f(1)*padic_(7,-2,12,-18446744073709551617) + *ORDEROLD*padic_(7,3,12,4294967297)*g(2); +Local ORDERG = g(2)*padic_(7,3,12,4294967297) + *padic_(7,0,12,4294967295)*f(1)*padic_(7,-2,12,-18446744073709551617); +Local ORDERZERO = ORDERF-ORDERG; +Local ORDERINVALID = padic_(7,0,12,3/2)+padic_(7,0,12,4294967295,1); +.sort +.sort +Print; +.end +#require flint? +#pend_if wordsize == 2 +assert succeeded? +assert result("STOREDCOMPACT") =~ expr("padic_(7,0,12,4294967295)") +assert result("STOREDPOSITIVE") =~ expr("padic_(7,0,12,4294967295)") +assert result("STOREDNEGATIVE") =~ expr("padic_(7,0,12,-4294967295)") +assert result("STOREDNESTED") =~ expr("f(padic_(7,0,12,18446744073709551617))") +assert result("STOREDSMALL") =~ expr("padic_(7,0,12,-2147483648)") + +assert result("ORDERZERO") =~ expr("0") +assert result("ORDERF") =~ expr(result("ORDERG")) +assert result("ORDERF") =~ expr("padic_(7,-2,12,-18446744073709551617) + *padic_(7,0,12,4294967295)*padic_(7,3,12,4294967297)*f(1)*g(2)") +assert result("ORDERINVALID") =~ expr("padic_(7,0,12,3/2)+padic_(7,0,12,4294967295,1)") +*--#] padic_inactive_algebra : +*--#[ padic_dummy_indices : +Indices i,j; +CFunctions f,g; +#StartPadic 7,100 +Global UNIT = padic_(7,0,100,340282366920938463463374607431768211455); +.sort +#EndPadic +Local F = UNIT*(f(i)*g(i)-f(j)*g(j)); +Sum i,j; +Print; +.end +#require flint? +assert succeeded? +assert result("F") =~ expr("0") +*--#] padic_dummy_indices : +*--#[ padic_save_load : +Symbols x; +CFunctions f; +#StartPadic 7,100 +Global UNIT = padic_(7,0,100,340282366920938463463374607431768211455); +.sort +#EndPadic +Global F = x*UNIT + + f(padic_(7,-2,100,-18446744073709551617)); +.store +Save compact.sav,F; +.end +Symbols y,x; +CFunctions g,f; +Load compact.sav; +Local G = F; +id x = 2; +Print; +.end +#require flint? +assert succeeded? +assert result("G") =~ expr("2*padic_(7,0,100,340282366920938463463374607431768211455) + + f(padic_(7,-2,100,-18446744073709551617))") +*--#] padic_save_load : +*--#[ padic_save_load_portability : +* Convert compact units across WORD sizes and endianness. +* Construct Save files with no user variables. Two units exercise width +* conversion followed by an enclosing tail, and a separate endian conversion. +Format 255; +Load compact16le.sav; +Local WIDTHPOSITIVE = F; +Print; +.end +Format 255; +Load compact32be.sav; +Local WIDTHNEGATIVE = F; +Print; +.end + +* Loading ordinary units must not expand their representation. +* Eleven ordinary units fit in 169 WORDs; inactive normalization must not +* grow them into compact wrappers requiring 213 WORDs. +#: maxtermsize 200 +Load compact-growth.sav; +Local GROWG = F; +Print; +.end +#require flint? +#pend_if wordsize == 2 +# Save/load conversion writes uninitialized metadata and index padding. +# Remove this condition once the standalone save/load fix is merged. +#pend_if valgrind? +#prepare [["16le", "<", 1, 2], ["32be", ">", -1, 4]].each do |suffix, endian, sign, width| +#prepare word = "#{width == 2 ? 's' : 'l'}#{endian}"; position = 1024+4*width +#prepare header = ([255]*8+[width,8,8,width == 2 ? 8 : 4]).pack("C*") +#prepare header += (endian == "<" ? [4,3,2,1] : [1,2,3,4]).pack("C*").ljust(16, "\0") +#prepare header += [48,40,40,72].pack("C*") +#prepare header += [10000].pack(word).ljust(16, "\0")+[6100].pack(word).ljust(16, "\0") +#prepare header += [2].pack("C"); header = header.ljust(512, "\0") +#prepare units = width == 2 ? [[-1]*8, [1,0,0,0,1]] : [[-1]*4, [1,0,1]] +#prepare functions = units.each_with_index.map do |limbs, i| +#prepare count = limbs.length; unit = [count+9,0,count+7,13,count+3,sign*(i == 0 ? count : -count)]+limbs+[1,1,3] +#prepare [127,9+unit.length,0,-16,7,-16,0,-16,100]+unit +#prepare end +#prepare body = functions.flatten; term = [body.length+4]+body+[1,1,3] +#prepare data = ([9,6,5,0,1,0,1,1,3]+term+[0]).pack("#{word}*") +#prepare entry = [position,data.bytesize,1024,0].pack("q#{endian}*")+([0]*5).pack("#{word}*") +#prepare entry += "F".ljust(17, "\0"); entry = entry.ljust(width == 2 ? 64 : 72, "\0") +#prepare index = ([0,1].pack("q#{endian}*")+entry).ljust(512, "\0") +#prepare write "compact#{suffix}.sav", header+index+"\0"*(4*width)+data +#prepare end +#prepare endian = [1].pack("l").bytes.first == 1 ? "<" : ">"; word = "s#{endian}" +#prepare header = ([255]*8+[2,8,8,8]).pack("C*") +#prepare header += (endian == "<" ? [4,3,2,1] : [1,2,3,4]).pack("C*").ljust(16, "\0") +#prepare header += [48,40,40,72].pack("C*") +#prepare header += [10000].pack(word).ljust(16, "\0")+[6100].pack(word).ljust(16, "\0") +#prepare header += [2].pack("C"); header = header.ljust(512, "\0") +#prepare body = (0..10).flat_map { |v| [127,17,0,-16,7,-16,v,-16,12,8,0,6,-1,-1,1,0,5] } +#prepare term = [body.length+4]+body+[1,1,3] +#prepare data = ([9,6,5,0,1,0,1,1,3]+term+[0]).pack("#{word}*") +#prepare entry = [1032,data.bytesize,1024,0].pack("q#{endian}*")+([0]*5).pack("#{word}*") +#prepare entry += "F".ljust(17, "\0"); entry = entry.ljust(64, "\0") +#prepare index = ([0,1].pack("q#{endian}*")+entry).ljust(512, "\0") +#prepare write "compact-growth.sav", header+index+"\0"*8+data +assert succeeded? +assert result("WIDTHPOSITIVE") =~ expr("padic_(7,0,100,-18446744073709551617)*padic_(7,0,100,340282366920938463463374607431768211455)") +assert result("WIDTHNEGATIVE") =~ expr("padic_(7,0,100,18446744073709551617)*padic_(7,0,100,-340282366920938463463374607431768211455)") + +expected = (0..10).map { |v| "padic_(7,#{v},12,4294967295)" }.join("*") +assert result("GROWG") =~ expr(expected) +*--#] padic_save_load_portability : +*--#[ padic_checkpoint_inactive : +#system -e test -f compact-start.tmp +#define first "`SYSTEMERROR_'" +#if `first' + #write "" + #close + #terminate 0 +#endif +Symbols x; +CFunctions f; +#StartPadic 7,100 +Global UNIT = padic_(7,0,100,340282366920938463463374607431768211455); +.sort +#EndPadic +On checkpoint; +Local F = x*UNIT + + f(padic_(7,-2,100,-18446744073709551617)); +#system sleep 1 +.sort +#system -e test -f compact-resume.tmp +#define stop "`SYSTEMERROR_'" +#write "" +#close +#if `stop' + #terminate +#endif +#remove +id x = 2; +Print; +.end +#system -e `FORM' 1.frm +#system `FORM' -R 1.frm +.end +#require flint? +#require unix? +#pend_if mpi? || valgrind? +assert return_value == 0 +assert result("F") =~ expr("2*padic_(7,0,100,340282366920938463463374607431768211455) + + f(padic_(7,-2,100,-18446744073709551617))") +*--#] padic_checkpoint_inactive : +*--#[ padic_storage_size : +#StartPadic 7,2000 +Format PadicPrint off; +CFunctions f; +Global F = padic_(7,0,2000,2^4096-1); +Global COMPACT = padic_(7,0,12,4294967295); +Print; +.sort +#EndPadic +Global ORDINARY = padic_(7,0,2000,2^4096-1); +Global POSITIVE = padic_(7,0,12,4294967295); +Global NEGATIVE = padic_(7,0,12,-4294967295); +Global NESTED = f(padic_(7,0,12,18446744073709551617)); +Global SMALL = padic_(7,0,12,-2147483648); +Print; +.store +Save integer.sav,POSITIVE,NEGATIVE,NESTED,SMALL; +Save compact.sav,COMPACT; +.end +#require flint? +assert succeeded? +large = expr("padic_(7,0,2000,#{2**4096-1})") +2.times do |pass| + assert result("F", pass).gsub(/\\\s*/, "") =~ large +end +assert result("ORDINARY").gsub(/\\\s*/, "") =~ large +# Statistics include the n+22 WORD bare term and a terminating zero WORD. +limbs = (4096+8*wordsize-1)/(8*wordsize) +2.times do |pass| + assert bytesize("F", pass) == (limbs+23)*wordsize +end +assert bytesize("ORDINARY") == (2*limbs+18)*wordsize +# Check the private Save encoding only on its original WORD size. +# Inactive values, including nested and short units, keep ordinary arguments. +if wordsize == 4 + integer = File.binread(File.join(workdir, "integer.sav")) + assert integer.include?([4,-1,1,3].pack("l*")) + assert integer.include?([4,-1,1,-3].pack("l*")) + assert integer.include?([8,1,0,1,1,0,0,7].pack("l*")) + assert !integer.include?([13,4,1,-1,1,1,3].pack("l*")) + assert !integer.include?([13,4,-1,-1,1,1,3].pack("l*")) + assert !integer.include?([13,6,3,1,0,1,1,1,3].pack("l*")) + assert integer.include?([-16,-2147483648].pack("l*")) + compact = File.binread(File.join(workdir, "compact.sav")) + assert compact.include?([13,4,1,-1,1,1,3].pack("l*")) +end +*--#] padic_storage_size : +*--#[ padic_normalization_size_limit : +* A normalized term exactly at MaxTermSize fits. +* The packed coefficient and f's arguments occupy exactly 200 WORDs. +#: maxtermsize 200 +#StartPadic 7,1080 +CFunctions f; +Format 255; +Format PadicPrint off; +Local FITF = f(<1>,...,<40>)*padic_(7,0,1080,-1); +Print; +.end +#require flint? +#pend_if wordsize == 2 +assert succeeded? +output = result("FITF").gsub(/\\\s*/, "") +assert output =~ expr("f(#{(1..40).to_a.join(",")})*padic_(7,0,1080,#{7**1080-1})") +*--#] padic_normalization_size_limit : +*--#[ padic_normalization_size_error : +* One more argument exceeds MaxTermSize. +#: maxtermsize 200 +#StartPadic 7,1080 +CFunctions f; +Local OVERF = f(<1>,...,<41>)*padic_(7,0,1080,-1); +.end +#require flint? +#pend_if wordsize == 2 +assert runtime_error?("Term too complex during normalization") +assert stdout.scan("Term too complex during normalization").size == 1 +*--#] padic_normalization_size_error : +*--#[ padic_scaling_size_error : +#: maxtermsize 200 +#StartPadic 7,1080 +CFunctions f; +Local F = -2*f(<1>,...,<41>)*padic_(7,0,1080,1); +.end +#require flint? +#pend_if wordsize == 2 +assert runtime_error?("Term too complex during normalization") +assert stdout.scan("Term too complex during normalization").size == 1 +*--#] padic_scaling_size_error : +*--#[ padic_normalization_scratch_error : +#: maxtermsize 200 +#StartPadic 7,2200 +Local F = padic_(7,0,2200,-1); +.end +#require flint? +#pend_if wordsize == 2 +assert runtime_error?("Term too complex during normalization") +assert stdout.scan("Term too complex during normalization").size == 1 +*--#] padic_normalization_scratch_error : +*--#[ padic_from_padic_size_error : +#: maxtermsize 200 +#StartPadic 7,1090 +CFunctions f; +Local F = -1; +ToPadic; +.sort +FromPadic,f; +.end +#require flint? +#pend_if wordsize == 2 +assert runtime_error?("Term too complex during normalization") +assert stdout !~ /Term too complex after FromPadic\./ +*--#] padic_from_padic_size_error : +*--#[ padic_prime_mismatch : +#StartPadic 7,10 +Format PadicPrint off; +Local F = padic_(11,0,10,1); +.end +#require flint? +#pend_if wordsize == 2 +assert runtime_error?("Prime in padic_ does not match the active p-adic prime.") +assert return_value != 0 +*--#] padic_prime_mismatch : +*--#[ padic_prime_mismatch_zero : +#StartPadic 7,10 +Local F = padic_(11,0,10,0); +.end +#require flint? +#pend_if wordsize == 2 +assert runtime_error?("Prime in padic_ does not match the active p-adic prime.") +assert return_value != 0 +*--#] padic_prime_mismatch_zero : +*--#[ padic_prime_mismatch_large : +#StartPadic 170141183460469231731687303715884105727,2 +Format PadicPrint off; +Local F = padic_(618970019642690137449562111,0,2,1); +.end +#require flint? +#pend_if wordsize == 2 +assert runtime_error?("Prime in padic_ does not match the active p-adic prime.") +assert return_value != 0 +*--#] padic_prime_mismatch_large : +*--#[ padic_prime_switch_mismatch : +#StartPadic 7,10 +Format PadicPrint off; +Local F = 31; +ToPadic; +.sort +#EndPadic +#StartPadic 11,10 +Format PadicPrint off; +Local G = F; +.end +#require flint? +#pend_if wordsize == 2 +assert runtime_error?("Prime in padic_ does not match the active p-adic prime.") +assert return_value != 0 +*--#] padic_prime_switch_mismatch : +*--#[ padic_inert_without_flint : +Global F = padic_(7,0,1000,2^4096-1); +Print; +.store +Save inert.sav,F; +.end +Load inert.sav; +Local G = F; +Print; +.end +#require !flint? +assert succeeded? +large = expr("padic_(7,0,1000,#{2**4096-1})") +assert result("F").gsub(/\\\s*/, "") =~ large +assert result("G").gsub(/\\\s*/, "") =~ large +# Without FLINT the literal keeps the ordinary padded integer argument. +limbs = (4096+8*wordsize-1)/(8*wordsize) +assert bytesize("F") == (2*limbs+18)*wordsize +*--#] padic_inert_without_flint : *--#[ humanstats : #- On humanstats; diff --git a/configure.ac b/configure.ac index ae64ba56e..364f8edfe 100644 --- a/configure.ac +++ b/configure.ac @@ -393,7 +393,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" @@ -423,7 +423,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 d8c3caeff..4bda768e2 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 12a2308fe..e694504e0 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 4c281613b..07451f1c4 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 ) ) { @@ -2069,6 +2105,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; } @@ -2078,6 +2127,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; } @@ -3604,6 +3659,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 == '-' ) { @@ -5403,6 +5464,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); } @@ -5435,6 +5502,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 61e45c55a..0cc0b75a5 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 cd4fb2e40..7a68458a2 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 : /* @@ -2799,6 +2852,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; } @@ -4123,6 +4204,33 @@ 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 ) { + if ( withpadic == 1 ) { + padicret = NormalizePadicUnit(BHEAD padicaccum,padicfunction); + if ( padicret == 1 ) goto NormZero; + if ( padicret > 1 ) padicfunction = padicaccum; + } + 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_ : @@ -4182,6 +4290,9 @@ regularratfun:; AT.NormDepth--; TermFree(n_llnum,"n_llnum"); TermFree(n_coef,"NormCoef"); +#ifdef WITHFLINT + TermFree(padicaccum,"NormPadic"); +#endif return(1); } else { @@ -4210,6 +4321,9 @@ regularratfun:; AT.NormDepth--; TermFree(n_llnum,"n_llnum"); TermFree(n_coef,"NormCoef"); +#ifdef WITHFLINT + TermFree(padicaccum,"NormPadic"); +#endif return(regval); NormInf: @@ -4236,12 +4350,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: @@ -4251,6 +4371,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 000000000..6e6d66eb8 --- /dev/null +++ b/sources/padic.c @@ -0,0 +1,2006 @@ +/** @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 +#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. + + The internal `padic_` coefficient carrier stores p explicitly. + Its prime must match the active `PadicContext` before the coefficient + is used in p-adic arithmetic. +*/ +#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; + /* Scratch rational used by FORM <-> padic/fmpq conversions. */ + fmpq_t q1; + /* Scratch integer used to check the prime stored in padic_. */ + fmpz_t prime; +} 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 pauxq1 (PadicAux->q1) +#define pauxprime (PadicAux->prime) + +static WORD GetPadicInteger(fmpz_t in, WORD *out); +static void SetPadicInteger(fmpz_t out, UWORD *in, WORD count); +static int TestNativePadicShape(const WORD *fun); +static int UnpackedPadicKnownReduced(const padic_t in, const WORD *arg); +static void SetPadicRational(padic_t out, const fmpq_t in); +static int SmallPadicRationalFits(slong N, const fmpq_t q); +static void DivideSmallPadicRational(padic_t out, const fmpq_t q, int bounded); +static const fmpz *CachedPadicModulus(PHEAD slong N, slong v); +static void AddReducedPadics(PHEAD padic_t out, const padic_t in1, const padic_t in2); +static int AddPadics(PHEAD WORD *outfun, WORD *in1, WORD *in2); +static int AddRatToPadic(PHEAD WORD *outfun, WORD *infun, UWORD *formrat, WORD nrat); + +/* + #] 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); + fmpq_init(aux->q1); + fmpz_init(aux->prime); +} +/* + #] InitPadicAux : + #[ ClearSinglePadicAux : +*/ +static void ClearSinglePadicAux(PADIC_AUX *aux) +{ + padic_clear(aux->p1); + padic_clear(aux->p2); + padic_clear(aux->p3); + fmpq_clear(aux->q1); + fmpz_clear(aux->prime); +} +/* + #] ClearSinglePadicAux : + #[ AllocatePadicAux : +*/ +static void AllocatePadicAux(void) +{ +#ifdef WITHPTHREADS + int id, totnum = AM.totalnumberofthreads; + totnum += GetNumberOfSortBots(totnum - 1); + 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 = AM.totalnumberofthreads; + /* Startup errors can precede allocation of the thread-private structures. */ + if ( AB == NULL || AB[0] == NULL ) return; + totnum += GetNumberOfSortBots(totnum - 1); + 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. + Include p^N when representable. For N < 0 an empty range is sufficient + and FLINT computes any needed positive powers on demand. + */ + ActivePadicContext = Malloc1(sizeof(padic_ctx_struct),"PadicContext"); + ctxmax = (N > 0) ? (slong)(N+(N < (FLINT_BITS == 64 ? INT64_MAX : INT32_MAX))) : 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 : + #[ GetPadicInteger : + + Converts a FLINT integer to the internal FORM integer encoding. + On the supported limb layout, copy positive integers directly from a + read-only GMP view. Small integers, negative integers and other layouts + use the shared conversion, including its established sign handling. +*/ +static WORD GetPadicInteger(fmpz_t in, WORD *out) +{ +#if BITSINWORD == 32 && FLINT_BITS == 64 && GMP_NUMB_BITS == 64 && GMP_NAIL_BITS == 0 && defined(mpz_limbs_read) + slong n = fmpz_size(in), i; + if ( n > 2 && fmpz_sgn(in) > 0 ) { + mpz_t view; + mp_srcptr limbs; + WORD count = (WORD)(2*n); + flint_mpz_init_set_readonly(view,in); + limbs = mpz_limbs_read(view); + for ( i = 0; i < n; i++ ) { + out[2*i] = (WORD)(UWORD)limbs[i]; + out[2*i+1] = (WORD)(UWORD)(limbs[i] >> BITSINWORD); + } + flint_mpz_clear_readonly(view); + if ( out[count-1] == 0 ) count--; + return(count); + } +#endif + return(flint_fmpz_get_form(in,out)); +} +/* + #] GetPadicInteger : + #[ SetPadicInteger : + + Converts the internal FORM integer encoding to a FLINT integer. + On the supported limb layout, combine at most 256 FORM words in a fixed + 128-limb (1 KiB) buffer. Other layouts and sizes use the shared conversion. +*/ +static void SetPadicInteger(fmpz_t out, UWORD *in, WORD count) +{ +#if BITSINWORD == 32 && FLINT_BITS == 64 && GMP_NUMB_BITS == 64 && GMP_NAIL_BITS == 0 + LONG n = count, i; + if ( n < 0 ) n = -n; + if ( n > 4 && n <= 256 ) { + ulong limbs[128]; + while ( n > 0 && in[n-1] == 0 ) n--; + if ( n > 4 ) { + for ( i = 0; i < (n+1)/2; i++ ) { + limbs[i] = (ulong)in[2*i]; + if ( 2*i+1 < n ) limbs[i] |= (ulong)in[2*i+1] << BITSINWORD; + } + fmpz_set_ui_array(out,limbs,(slong)((n+1)/2)); + if ( count < 0 ) fmpz_neg(out,out); + return; + } + } +#endif + flint_fmpz_set_form(out,in,count); +} +/* + #] SetPadicInteger : + #] 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_(p, v, N, u) + where + - p is the prime, + - v is the p-adic valuation, + - N is the precision stored with this coefficient, + - u is the FLINT unit. + + The prime p must match the active p-adic context configured through + #StartPadic / #EndPadic. Both p and u are stored as FORM integers with + denominator 1 on input. Inactive normalization preserves the unit encoding. + Active arithmetic packs units as -SNUMBER or a private PADICLIMBS object + in a one-term argument with coefficient 1/1. Its positive length + is PADICLIMBHEAD+n; its signed count is +/-n, followed by n unsigned + FORM limbs, least significant first. Zero uses -SNUMBER. Recognition and + FromPadic do not require an active context. + + #] Explanations : + #[ TestPadicLimbTerm : +*/ +int TestPadicLimbTerm(const WORD *term) +{ + LONG n; + if ( term[0] < PADICLIMBHEAD+5 || term[1] != PADICLIMBS ) return(0); + if ( term[2] != term[0]-4 || term[3] == WORD_MIN_VALUE ) return(0); + n = ABS(term[3]); + if ( n == 0 || term[2] != PADICLIMBHEAD+n ) return(0); + if ( (UWORD)term[PADICLIMBHEAD+n] == 0 ) return(0); + return(term[term[0]-3] == 1 && term[term[0]-2] == 1 + && term[term[0]-1] == 3); +} +/* + #] TestPadicLimbTerm : + #[ ExpandPadicArgument : + Write a validated PADICLIMBS argument as an ordinary integer. + Callers must provide enough output space and ensure the length fits in WORD. +*/ +int ExpandPadicArgument(WORD *out, const WORD *arg) +{ + const WORD *limbs = arg+ARGHEAD+1+PADICLIMBHEAD; + WORD count = arg[ARGHEAD+PADICLIMBHEAD]; + LONG n = ABS(count), size = ARGHEAD+2*n+2, i; + *out++ = (WORD)size; + *out++ = 0; + FILLARG(out); + *out++ = (WORD)(2*n+2); + i = n; NCOPY(out,limbs,i); + *out++ = 1; + for ( i = 1; i < n; i++ ) *out++ = 0; + *out = (WORD)(count < 0 ? -(2*n+1) : 2*n+1); + return((int)size); +} +/* + #] ExpandPadicArgument : + #[ TestNativePadicShape : + + Checks the common padic_(p,v,N,u) encoding with WORD-sized p, v and N. + The unit may be a small number or a compact limb term. Return 1 for a + valid shape, otherwise 0 so TestPadic performs its general checks. + This checks the encoding only; it does not certify a reduced unit. +*/ +static int TestNativePadicShape(const WORD *fun) +{ + const WORD *arg; + WORD length = fun[1]; + if ( length < FUNHEAD+8 ) return(0); + arg = fun+FUNHEAD; + if ( arg[0] != -SNUMBER || arg[1] <= 1 + || arg[2] != -SNUMBER || arg[4] != -SNUMBER ) return(0); + arg += 6; + if ( arg[0] == -SNUMBER ) return(length == FUNHEAD+8); + if ( length < FUNHEAD+6+ARGHEAD+PADICLIMBHEAD+5 + || arg[0] != length-FUNHEAD-6 + || arg[ARGHEAD] != arg[0]-ARGHEAD ) return(0); + return(TestPadicLimbTerm(arg+ARGHEAD)); +} +/* + #] TestNativePadicShape : + #[ TestPadic : + + Checks whether fun has the shape of a legal padic_ with its arguments: + padic_(p,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; + if ( TestNativePadicShape(fun) ) return(1); + f = fun + FUNHEAD; + fstop = fun + fun[1]; + nargs = 0; + while ( f < fstop ) { + nargs++; + NEXTARG(f); + } + if ( nargs != 4 || f != fstop ) return(0); + + f = fun + FUNHEAD; + /* p: integer greater than 1. Matching the active prime is checked separately. */ + if ( *f == -SNUMBER ) { + if ( f[1] <= 1 ) return(0); + f += 2; + } + else { + if ( *f < ARGHEAD+4 || f[*f-1] <= 0 ) return(0); + nnum = (WORD)((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); + if ( f[ARGHEAD+nnum+1] != 1 ) return(0); + for ( i = 1; i < nnum; i++ ) { + if ( f[ARGHEAD+nnum+1+i] != 0 ) return(0); + } + if ( (UWORD)f[ARGHEAD+1] <= 1 ) { + for ( i = 1; i < nnum; i++ ) { + if ( f[ARGHEAD+1+i] != 0 ) break; + } + if ( i == nnum ) return(0); + } + f += *f; + } + /* 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 or private wrapped limb-list. */ + if ( *f == -SNUMBER ) { + f += 2; + } + else if ( *f > 0 && f[ARGHEAD] == *f-ARGHEAD && TestPadicLimbTerm(f+ARGHEAD) ) { + return(1); + } + else { + nnum = (WORD)((ABS(fstop[-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 : + #[ CheckPadicPrime : + + Checks the prime of a structurally valid padic_ against the active context. + A mismatch is fatal. With no active context there is nothing to compare. +*/ +void CheckPadicPrime(WORD *fun) +{ + WORD *f = fun + FUNHEAD; + int matches; + + if ( !PadicActive ) return; + if ( *f == -SNUMBER ) { + matches = fmpz_equal_si(PadicContext->p,(slong)f[1]); + } + else { + WORD nnum = (WORD)((f[*f-1]-1)/2); + GETIDENTITY + + flint_fmpz_set_form(pauxprime,(UWORD *)(f+ARGHEAD+1),nnum); + matches = fmpz_equal(pauxprime,PadicContext->p); + } + if ( !matches ) { + MLOCK(ErrorMessageLock); + MesPrint("Prime in padic_ does not match the active p-adic prime."); + MUNLOCK(ErrorMessageLock); + Terminate(-1); + } +} +/* + #] CheckPadicPrime : + #[ UnpackPadic : + + Converts the arguments of a padic_ function back into FLINT's padic_t. + Callers must ensure that fun is a legal padic_ with decodable arguments, + either by a successful TestPadic() call or by using output from PackPadic(). + The p-adic system must be active. +*/ +static void 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); + } + + CheckPadicPrime(fun); + + /* + Skip the checked prime, then read v, N, u. + */ + f = fun + FUNHEAD; + NEXTARG(f); + + 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])); + } + /* An ordinary integer's term consists only of its coefficient. */ + else if ( f[ARGHEAD] == ABS(f[*f-1])+1 ) { + WORD size = f[*f-1], nnum = (WORD)((ABS(size)-1)/2); + while ( nnum && f[ARGHEAD+nnum] == 0 ) nnum--; + SetPadicInteger(padic_unit(out),(UWORD *)(f+ARGHEAD+1), + (size < 0) ? -nnum : nnum); + } + else { + SetPadicInteger(padic_unit(out),(UWORD *)(f+ARGHEAD+1+PADICLIMBHEAD), + f[ARGHEAD+PADICLIMBHEAD]); + } + if ( UnpackedPadicKnownReduced(out,f) ) return; + padic_reduce(out,PadicContext); +} +/* + #] UnpackPadic : + #[ PackPadic : + + Converts a FLINT's padic_t into FORM's internal padic notation: + padic_(p,v,N,u) + Callers must provide enough output space for the packed function. + The FLINT value must already be reduced: for nonzero values, p does not + divide u and 0 < u < p^(N-v); zero has u = v = 0. + + 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; + + v = (LONG)padic_val(in); + N = (LONG)padic_prec(in); + + t = fun; + *t++ = PADICFUN; + t++; // fun[1] (function length) is filled at the end. + FILLFUN(t); + + if ( fmpz_cmp_ui(PadicContext->p,WORD_MAX_VALUE) <= 0 ) { + *t++ = -SNUMBER; + *t++ = (WORD)fmpz_get_si(PadicContext->p); + } + else { + arg = t + ARGHEAD + 1; + nnum = GetPadicInteger(PadicContext->p,arg); + arg += nnum; + *arg++ = 1; + for ( i = 1; i < nnum; i++ ) *arg++ = 0; + *arg++ = 2*nnum+1; + *t++ = ARGHEAD + 2*nnum + 2; + *t++ = 0; + FILLARG(t); + *t++ = 2*nnum + 2; + t = arg; + } + + /* + Pack valuation v and N: + compact -SNUMBER when it fits in WORD, otherwise as a two-word + long integer argument. + */ + for ( i = 0; i < 2; i++ ) { + LONG value = i == 0 ? v : N; + if ( value >= WORD_MIN_VALUE && value <= WORD_MAX_VALUE ) { + *t++ = -SNUMBER; + *t++ = (WORD)value; + } + else { + x = LongAbs(value); + *t++ = ARGHEAD + 6; + *t++ = 0; + FILLARG(t); + *t++ = 6; + *t++ = (WORD)(UWORD)x; + *t++ = (WORD)(UWORD)(x >> BITSINWORD); + *t++ = 1; + *t++ = 0; + *t++ = value < 0 ? -5 : 5; + } + } + if ( fmpz_fits_si(padic_unit(in)) + && (small = (LONG)fmpz_get_si(padic_unit(in)), + small >= WORD_MIN_VALUE && small <= WORD_MAX_VALUE) ) { + *t++ = -SNUMBER; + *t++ = (WORD)small; + } + else { + arg = t + ARGHEAD + 1 + PADICLIMBHEAD; + nnum = GetPadicInteger(padic_unit(in),arg); + nabs = ABS(nnum); + t[0] = ARGHEAD + 1 + PADICLIMBHEAD + nabs + 3; + for ( i = 1; i < ARGHEAD; i++ ) t[i] = 0; + t[ARGHEAD] = t[0] - ARGHEAD; + t[ARGHEAD+1] = PADICLIMBS; + t[ARGHEAD+2] = PADICLIMBHEAD + nabs; + t[ARGHEAD+3] = nnum; + t += *t; + t[-3] = 1; + t[-2] = 1; + t[-1] = 3; + } + fun[1] = (WORD)(t - fun); + return(fun[1]); +} +/* + #] PackPadic : + #[ PadicUnitKnownReduced : + + Certify a positive unit using exact bounds on p^k and divisibility. + Near powers of two, p^k >= (2^b-k*(2^b-p))*2^(b*(k-1)). + Uncertain bounds retain FLINT; imported units may reuse GMP's remainder. +*/ +static int PadicUnitKnownReduced(const UWORD *limbs, WORD count, LONG v, LONG N, const fmpz *imported) +{ + ULONG p, k, bits, pbits, base, delta, lower, shift, top, rem = 0; + LONG i, offset; + int inrange; + if ( count <= 0 || v >= N || fmpz_cmp_ui(PadicContext->p,WORD_MAX_VALUE) > 0 ) + return(0); + p = (ULONG)fmpz_get_ui(PadicContext->p); + k = (ULONG)N-(ULONG)v; + bits = (ULONG)(count-1)*BITSINWORD+FLINT_BIT_COUNT(limbs[count-1]); + pbits = FLINT_BIT_COUNT(p)-1; + inrange = k >= bits || k*pbits >= bits; + if ( !inrange ) { + pbits = FLINT_BIT_COUNT(p-1); + base = (ULONG)1 << pbits; + delta = base-p; + if ( delta == 0 || k <= (base-1)/delta ) { + lower = base-k*delta; + shift = pbits*(k-1); + top = shift+FLINT_BIT_COUNT(lower); + if ( bits < top ) inrange = 1; + else if ( bits == top ) { + i = (LONG)(shift/BITSINWORD); offset = (LONG)(shift%BITSINWORD); + top = (ULONG)limbs[i] >> offset; + if ( offset && i+1 < count ) + top |= (ULONG)limbs[i+1] << (BITSINWORD-offset); + inrange = top < lower; + } + } + } + if ( !inrange ) return(0); + if ( BITSINWORD == 32 && FLINT_BITS == 64 && GMP_NUMB_BITS == 64 && GMP_NAIL_BITS == 0 + && imported && count >= 12 ) return(fmpz_fdiv_ui(imported,(ulong)p) != 0); + for ( i = count-1; i >= 0; i-- ) rem = ((rem << BITSINWORD)+limbs[i]) % p; + return(rem != 0); +} +/* + #] PadicUnitKnownReduced : + #[ UnpackedPadicKnownReduced : + + Checks whether the unit imported from a validated padic_ carrier is + already reduced. Reuse its FORM limbs for the range bound and its + FLINT integer for the divisibility check when useful. Return 1 only + when both checks succeed within FLINT's supported exponent range; + otherwise UnpackPadic must call padic_reduce. +*/ +static int UnpackedPadicKnownReduced(const padic_t in, const WORD *arg) +{ + const UWORD *limbs; + UWORD small; + WORD count; + slong N = padic_prec(in), v = padic_val(in); + if ( N <= v || (ulong)N-(ulong)v + > (ulong)(FLINT_BITS == 64 ? INT64_MAX : INT32_MAX) ) return(0); + if ( *arg == -SNUMBER ) { + if ( arg[1] <= 0 ) return(0); + small = (UWORD)arg[1]; limbs = &small; count = 1; + } + else if ( arg[ARGHEAD] == ABS(arg[*arg-1])+1 ) { + WORD size = arg[*arg-1]; + if ( size < 0 ) return(0); + count = (WORD)((size-1)/2); + limbs = (const UWORD *)(arg+ARGHEAD+1); + while ( count && limbs[count-1] == 0 ) count--; + } + else { + count = arg[ARGHEAD+PADICLIMBHEAD]; + limbs = (const UWORD *)(arg+ARGHEAD+1+PADICLIMBHEAD); + } + return(PadicUnitKnownReduced(limbs,count,v,N,padic_unit(in))); +} +/* + #] UnpackedPadicKnownReduced : + #[ NormalizePadicUnit : + + Normalize a lone coefficient only after the rational coefficient and any + other p-adics have been combined. Leave an unchanged compact unit in place. + Callers must have validated fun with TestPadic. + Return 0 for unchanged, 1 for zero, or the packed length. +*/ +int NormalizePadicUnit(PHEAD WORD *out, WORD *fun) +{ + GETBIDENTITY + const WORD *arg = fun+FUNHEAD, *varg = fun+FUNHEAD; + const UWORD *limbs = NULL; + UWORD small; + WORD count = 0; + ULONG value; + LONG v, N; + int same = 0; + NEXTARG(arg); NEXTARG(arg); NEXTARG(arg); + NEXTARG(varg); + if ( *varg == -SNUMBER ) v = varg[1]; + else { + value = ((ULONG)(UWORD)varg[ARGHEAD+2] << BITSINWORD) + + (UWORD)varg[ARGHEAD+1]; + v = ULongToLong(varg[ARGHEAD+5] < 0 ? -value : value); + } + NEXTARG(varg); + if ( *varg == -SNUMBER ) N = varg[1]; + else { + value = ((ULONG)(UWORD)varg[ARGHEAD+2] << BITSINWORD) + + (UWORD)varg[ARGHEAD+1]; + N = ULongToLong(varg[ARGHEAD+5] < 0 ? -value : value); + } + if ( *arg == -SNUMBER ) { + small = WordAbs(arg[1]); + limbs = &small; + count = small == 0 ? 0 : (arg[1] < 0 ? -1 : 1); + } + else if ( arg[ARGHEAD] != ABS(arg[*arg-1])+1 ) { + count = arg[ARGHEAD+PADICLIMBHEAD]; + limbs = (const UWORD *)(arg+ARGHEAD+1+PADICLIMBHEAD); + } + if ( count > 0 && PadicUnitKnownReduced(limbs,count,v,N,NULL) ) return(0); + UnpackPadic(paux1,fun); + if ( padic_is_zero(paux1) ) return(1); + if ( v == padic_val(paux1) && count > 0 ) { + if ( *arg == -SNUMBER ) same = fmpz_equal_si(padic_unit(paux1),arg[1]); +#if GMP_NUMB_BITS % BITSINWORD == 0 + else if ( (fmpz_bits(padic_unit(paux1))+BITSINWORD-1)/BITSINWORD == (ulong)count ) { + mpz_t unit; + LONG i; + flint_mpz_init_set_readonly(unit,padic_unit(paux1)); + same = 1; + for ( i = 0; i < count; i++ ) { + mp_limb_t limb = mpz_getlimbn(unit,(i*BITSINWORD)/GMP_NUMB_BITS); + if ( limbs[i] != (UWORD)(limb >> ((i*BITSINWORD)%GMP_NUMB_BITS)) ) { + same = 0; break; + } + } + flint_mpz_clear_readonly(unit); + } +#endif + } + if ( same ) return(0); + return(PackPadic(out,paux1)); +} +/* + #] NormalizePadicUnit : + #] 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 ) { + SetPadicInteger(fmpq_numref(result),num,sign*nnum); + } + else { + fmpz_zero(fmpq_numref(result)); + } + SetPadicInteger(fmpq_denref(result),den,nden); + + fmpq_canonicalise(result); +} +/* + #] FormRatToFmpq : + #[ SetPadicRational : + + Converts a canonical rational to a p-adic at the output precision. + Use the integer conversion or the small-fraction shortcut when possible. + Other fractions use the general FLINT conversion. +*/ +static void SetPadicRational(padic_t out, const fmpq_t in) +{ + if ( fmpz_is_one(fmpq_denref(in)) ) + padic_set_fmpz(out,fmpq_numref(in),PadicContext); + else if ( SmallPadicRationalFits(padic_prec(out),in) ) { + fmpz_set(padic_unit(out),fmpq_numref(in)); + padic_val(out) = 0; + DivideSmallPadicRational(out,in,1); + } + else padic_set_fmpq(out,in,PadicContext); +} +/* + #] SetPadicRational : + #[ SmallPadicRationalFits : + + Checks the precision and fraction bounds for the small-rational shortcut. +*/ +static int SmallPadicRationalFits(slong N, const fmpq_t q) +{ + return(N > 0 && !fmpz_is_zero(fmpq_numref(q)) + && fmpz_cmpabs(fmpq_numref(q),PadicContext->p) < 0 + && fmpz_cmp_ui(fmpq_denref(q),1) > 0 + && fmpz_cmp(fmpq_denref(q),PadicContext->p) < 0 + && fmpz_cmp_ui(fmpq_denref(q),WORD_MAX_VALUE) <= 0); +} +/* + #] SmallPadicRationalFits : + #[ DivideSmallPadicRational : + + Divides the prepared unit numerator by the small fraction's denominator + modulo p^(N-v). The caller has checked the fraction and set the valuation. +*/ +static void DivideSmallPadicRational(padic_t out, const fmpq_t q, int bounded) +{ + slong N = padic_prec(out), v = padic_val(out); + fmpz_t modulus; + ulong d, inverse, remainder, correction; + int owned; + d = fmpz_get_ui(fmpq_denref(q)); + owned = _padic_ctx_pow_ui(modulus,(ulong)(N-v),PadicContext); + /* Add a multiple of the modulus to make the numerator divisible by d. */ + inverse = n_invmod(fmpz_fdiv_ui(modulus,d),d); + remainder = fmpz_fdiv_ui(padic_unit(out),d); + correction = remainder ? d-n_mulmod2(remainder,inverse,d) : 0; + fmpz_addmul_ui(padic_unit(out),modulus,correction); + fmpz_divexact_ui(padic_unit(out),padic_unit(out),d); + /* A bounded proper fraction needs at most one modular correction. */ + if ( bounded && fmpz_cmpabs(fmpq_numref(q),fmpq_denref(q)) < 0 ) { + if ( fmpz_sgn(padic_unit(out)) < 0 ) + fmpz_add(padic_unit(out),padic_unit(out),modulus); + else if ( fmpz_cmp(padic_unit(out),modulus) >= 0 ) + fmpz_sub(padic_unit(out),padic_unit(out),modulus); + } + /* Otherwise use the general remainder operation. */ + else fmpz_mod(padic_unit(out),padic_unit(out),modulus); + if ( fmpz_is_zero(padic_unit(out)) ) padic_val(out) = 0; + if ( owned ) fmpz_clear(modulus); +} +/* + #] DivideSmallPadicRational : + #[ CachedPadicModulus : + + Returns the modulus p^(N-v) from the existing context cache, or from + the product of two cached powers when v is negative. Return NULL when + the required powers are unavailable. + + Uses the thread-local auxiliary integer as scratch. The returned pointer + is borrowed and must not be cleared or kept across another scratch use. +*/ +static const fmpz *CachedPadicModulus(PHEAD slong N, slong v) +{ + slong exponent; + if ( N <= v ) return(NULL); + /* This can potentially overflow. Should be fixed. */ + exponent = N-v; + if ( PadicContext->min <= exponent && exponent < PadicContext->max ) + return(PadicContext->pow+(exponent-PadicContext->min)); + if ( N <= 0 || v >= 0 || v < -N + || N < PadicContext->min || N >= PadicContext->max + || -v < PadicContext->min || -v >= PadicContext->max ) return(NULL); + fmpz_mul(pauxprime,PadicContext->pow+(N-PadicContext->min), + PadicContext->pow+(-v-PadicContext->min)); + return(pauxprime); +} +/* + #] CachedPadicModulus : + #[ AddReducedPadics : + + Adds two reduced p-adics at the output precision, using the cached + modulus when precision and valuation match. Other cases use FLINT. + Uses the thread-local auxiliary integer as scratch for the modulus. +*/ +static void AddReducedPadics(PHEAD padic_t out, const padic_t in1, const padic_t in2) +{ + slong N = FLINT_MIN(PadicPrecision,FLINT_MIN(padic_prec(in1),padic_prec(in2))); + slong v = padic_val(in1); + const fmpz *modulus; + padic_prec(out) = N; + /* The unit shortcut requires matching precision and valuation. */ + if ( v != padic_val(in2) || N != padic_prec(in1) || N != padic_prec(in2) ) goto fallback; + modulus = CachedPadicModulus(BHEAD N,v); + if ( !modulus ) goto fallback; + /* Each unit is below the modulus, so one subtraction suffices. */ + fmpz_add(padic_unit(out),padic_unit(in1),padic_unit(in2)); + if ( fmpz_cmp(padic_unit(out),modulus) >= 0 ) + fmpz_sub(padic_unit(out),padic_unit(out),modulus); + padic_val(out) = v; + /* Remove prime factors from the unit and update the valuation. */ + _padic_canonicalise(out,PadicContext); + return; +fallback: + padic_add(out,in1,in2,PadicContext); +} +/* + #] AddReducedPadics : + #[ AddRatToPadic : + + Adds a FORM rational coefficient to a padic_ function in infun. + The exact rational is converted at the input's precision, capped at the + active precision. Addition cannot recover discarded input digits. + + Return value: + - 0 on success, with a non-zero sum packed into outfun, + - 1 if the sum is p-adic zero, + - -1 on runtime/validation failure. +*/ +static int AddRatToPadic(PHEAD WORD *outfun, WORD *infun, UWORD *formrat, WORD nrat) +{ + if ( !PadicActive ) return(-1); + UnpackPadic(paux1,infun); + FormRatToFmpq(pauxq1,formrat,nrat); + padic_prec(paux2) = FLINT_MIN(PadicPrecision,padic_prec(paux1)); + SetPadicRational(paux2,pauxq1); + AddReducedPadics(BHEAD paux3,paux1,paux2); + PackPadic(outfun,paux3); + if ( padic_is_zero(paux3) ) return(1); + return(0); +} +/* + #] AddRatToPadic : + #[ 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) +{ + slong v, w, N; + if ( !PadicActive ) return(-1); + UnpackPadic(paux1,infun); + FormRatToFmpq(pauxq1,formrat,nrat); + v = padic_val(paux1); + N = FLINT_MIN(PadicPrecision,padic_prec(paux1)); + padic_prec(paux3) = N; + if ( padic_is_zero(paux1) || fmpq_is_zero(pauxq1) ) padic_zero(paux3); + else if ( fmpz_is_one(fmpq_denref(pauxq1)) && fmpz_equal_si(fmpq_numref(pauxq1),-1) ) { + padic_neg(paux3,paux1,PadicContext); + } + else if ( SmallPadicRationalFits(N,pauxq1) ) { + /* These exact fractions have valuation zero. Divide the unit at + N-v for either sign of v, retaining all required guard digits. */ + if ( v >= N ) padic_zero(paux3); + else { + fmpz_mul(padic_unit(paux3),padic_unit(paux1),fmpq_numref(pauxq1)); + padic_val(paux3) = v; + DivideSmallPadicRational(paux3,pauxq1,padic_prec(paux1) <= N); + } + } + else { + const fmpz *modulus; + /* An exact rational shifts the input error and valuation by w. */ + w = fmpz_remove(padic_unit(paux2),fmpq_numref(pauxq1),PadicContext->p); + w -= fmpz_remove(pauxprime,fmpq_denref(pauxq1),PadicContext->p); + /* This can potentially overflow. Should be fixed. */ + N = FLINT_MIN(PadicPrecision,padic_prec(paux1)+w); + v += w; + padic_prec(paux3) = N; + if ( v >= N ) padic_zero(paux3); + else { + /* Use the input valuation to retain the rational's guard digits. */ + padic_prec(paux2) = N-padic_val(paux1); + SetPadicRational(paux2,pauxq1); + /* Reduce the unit product using the cached modulus. */ + modulus = CachedPadicModulus(BHEAD N,v); + if ( modulus ) { + fmpz_mul(padic_unit(paux3),padic_unit(paux1),padic_unit(paux2)); + fmpz_mod(padic_unit(paux3),padic_unit(paux3),modulus); + padic_val(paux3) = padic_is_zero(paux3) ? 0 : v; + } + /* Use FLINT when the required modulus is not cached. */ + else padic_mul(paux3,paux1,paux2,PadicContext); + } + } + PackPadic(outfun,paux3); + if ( padic_is_zero(paux3) ) return(1); + return(0); +} +/* + #] MulRatToPadic : + #[ AddPadics : + + Adds two padic_ functions (in1 + in2) and stores the sum in outfun. + + Return value: + - 0 on success, with a non-zero sum packed into outfun, + - 1 if the sum is p-adic zero, + - -1 on runtime/validation failure. +*/ +static int AddPadics(PHEAD WORD *outfun, WORD *in1, WORD *in2) +{ + if ( !PadicActive ) return(-1); + UnpackPadic(paux1,in1); + UnpackPadic(paux2,in2); + AddReducedPadics(BHEAD paux3,paux1,paux2); + PackPadic(outfun,paux3); + if ( padic_is_zero(paux3) ) return(1); + return(0); +} +/* + #] AddPadics : + #[ 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) +{ + slong N1, N2; + if ( !PadicActive ) return(-1); + UnpackPadic(paux1,fun1); + UnpackPadic(paux2,fun2); + /* Preserve the smaller relative precision, shifted to the product. + This can potentially overflow. Should be fixed. */ + N1 = padic_prec(paux1)+padic_val(paux2); + N2 = padic_prec(paux2)+padic_val(paux1); + padic_prec(paux3) = FLINT_MIN(PadicPrecision,FLINT_MIN(N1,N2)); + 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) +{ + slong N1, N2, v; + if ( !PadicActive ) return(-1); + UnpackPadic(paux1,fun1); + UnpackPadic(paux2,fun2); + if ( padic_is_zero(paux2) ) { + MLOCK(ErrorMessageLock); + MesPrint("Division by zero in p-adic arithmetic."); + MUNLOCK(ErrorMessageLock); + Terminate(-1); + return(-1); + } + v = padic_val(paux2); + /* This can potentially overflow. Should be fixed. */ + N1 = padic_prec(paux1)-v; + N2 = padic_prec(paux2)+padic_val(paux1)-v-v; + padic_prec(paux3) = FLINT_MIN(PadicPrecision,FLINT_MIN(N1,N2)); + 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); + UnpackPadic(paux1,fun); + if ( padic_is_zero(paux1) ) { + MLOCK(ErrorMessageLock); + MesPrint("Division by zero in p-adic arithmetic."); + MUNLOCK(ErrorMessageLock); + Terminate(-1); + return(-1); + } + /* This can potentially overflow. Should be fixed. */ + padic_prec(paux3) = FLINT_MIN(PadicPrecision, + padic_prec(paux1)-padic_val(paux1)-padic_val(paux1)); + 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(p?,v?,N?,u?) = makerational_(u,p^(N-v))*p^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); + CheckPadicPrime(fun); + + // Skip the checked prime, then read v, N, u. + arg = fun + FUNHEAD; + NEXTARG(arg); + 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] = WordAbs(arg[1]); + } + else if ( arg[ARGHEAD] == ABS(arg[*arg-1])+1 ) { + nu = (WORD)((ABS(arg[*arg-1])-1)/2); + while ( nu && arg[ARGHEAD+nu] == 0 ) nu--; + if ( arg[*arg-1] < 0 ) nu = -nu; + for ( i = 0; i < ABS(nu); i++ ) u[i] = (UWORD)arg[ARGHEAD+1+i]; + } + else { + nu = arg[ARGHEAD+PADICLIMBHEAD]; + for ( i = 0; i < ABS(nu); i++ ) u[i] = (UWORD)arg[ARGHEAD+1+PADICLIMBHEAD+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 ( TestPadic(fun) == 0 ) return(0); + UnpackPadic(paux1,fun); + + 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) ) { + CheckPadicPrime(t); + return(Generator(BHEAD term,level)); + } + t += t[1]; /* advance to the next function in the term */ + } + } + + FormRatToFmpq(pauxq1,(UWORD *)tstop,ncoef); + padic_prec(paux1) = PadicPrecision; + SetPadicRational(paux1,pauxq1); + if ( padic_is_zero(paux1) ) return(0); + // Overwrite the old rational coefficient of the term by padic_(p,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_(p,v,N,u) at the current execution level is + converted to f(p,v,N,u), independently of the active context. +*/ +int FromPadic(PHEAD WORD *term, WORD level, WORD target) +{ + GETBIDENTITY + WORD *t, *first, *stop, *newterm, *out, *oldworkpointer = AT.WorkPointer; + int len, retval; + + stop = term+*term; + stop -= ABS(stop[-1]); + for ( t = term+1; t < stop; t += t[1] ) { + if ( *t == PADICFUN && TestPadic(t) ) break; + } + if ( t == stop ) return(Generator(BHEAD term,level)); + first = t; + if ( AT.WorkPointer < term+*term ) AT.WorkPointer = term+*term; + newterm = AT.WorkPointer; out = newterm+1; + t = term+1; len = first-t; + NCOPY(out,t,len); + while ( t < stop ) { + if ( t == first || (*t == PADICFUN && TestPadic(t)) ) { + WORD *arg = t+FUNHEAD, *fun = out; + NEXTARG(arg); NEXTARG(arg); NEXTARG(arg); + len = arg-t; NCOPY(out,t,len); + if ( *t == -SNUMBER ) { + len = 2; NCOPY(out,t,len); + } + else if ( t[ARGHEAD] == ABS(t[*t-1])+1 ) { + len = *t; NCOPY(out,t,len); + } + else { + len = ExpandPadicArgument(out,t); + out += len; + NEXTARG(t); + } + fun[1] = (WORD)(out-fun); + *fun = target; + } + else { + len = t[1]; NCOPY(out,t,len); + } + } + len = term+*term-stop; NCOPY(out,stop,len); + *newterm = (WORD)(out-newterm); + AT.WorkPointer = out; + retval = Generator(BHEAD newterm,level); + AT.WorkPointer = oldworkpointer; + return(retval); +} +/* + #] 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; + int padicret; + + 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); + + /* Compare1 validated the final object as the p-adic coefficient. */ + if ( AT.SortPadicMode == 3 ) { + /* Both coefficients are normalized padic_ records. */ + fun1 = s1+1; while ( fun1+fun1[1] < coef1 ) fun1 += fun1[1]; + fun2 = s2+1; while ( fun2+fun2[1] < coef2 ) fun2 += fun2[1]; + } + else if ( AT.SortPadicMode == 1 ) { + /* First coefficient is padic_, second is rational. */ + fun1 = s1+1; while ( fun1+fun1[1] < coef1 ) fun1 += fun1[1]; + fun2 = coef2; + } + else if ( AT.SortPadicMode == 2 ) { + /* Second coefficient is padic_, first is rational. */ + fun2 = s2+1; while ( fun2+fun2[1] < coef2 ) fun2 += fun2[1]; + fun1 = coef1; + } + else { + MLOCK(ErrorMessageLock); + MesPrint("Illegal value %d for AT.SortPadicMode in AddWithPadic.",AT.SortPadicMode); + MUNLOCK(ErrorMessageLock); + Terminate(-1); + return(0); + } + + fun3 = TermMalloc("AddWithPadic"); + if ( AT.SortPadicMode == 3 ) padicret = AddPadics(BHEAD fun3,fun1,fun2); + else if ( AT.SortPadicMode == 1 ) + padicret = AddRatToPadic(BHEAD fun3,fun1,(UWORD *)coef2,size2); + else padicret = AddRatToPadic(BHEAD fun3,fun2,(UWORD *)coef1,size1); + if ( padicret != 0 ) { + TermFree(fun3,"AddWithPadic"); + /* Terms cancel. A failed conversion leaves the inputs intact. */ + if ( padicret > 0 ) *ps1 = *ps2 = 0; + AT.SortPadicMode = 0; + return(0); + } + + 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, padicret; + + if ( !PadicActive ) return(0); + + coef1 = term1+*term1; size1 = coef1[-1]; coef1 -= ABS(size1); + coef2 = term2+*term2; size2 = coef2[-1]; coef2 -= ABS(size2); + /* Compare1 validated the final object as the p-adic coefficient. */ + if ( AT.SortPadicMode == 3 ) { + fun1 = term1+1; while ( fun1+fun1[1] < coef1 ) fun1 += fun1[1]; + fun2 = term2+1; while ( fun2+fun2[1] < coef2 ) fun2 += fun2[1]; + } + else if ( AT.SortPadicMode == 1 ) { + fun1 = term1+1; while ( fun1+fun1[1] < coef1 ) fun1 += fun1[1]; + fun2 = coef2; + } + else if ( AT.SortPadicMode == 2 ) { + fun2 = term2+1; while ( fun2+fun2[1] < coef2 ) fun2 += fun2[1]; + fun1 = coef1; + } + else { + MLOCK(ErrorMessageLock); + MesPrint("Illegal value %d for AT.SortPadicMode in MergeWithPadic.",AT.SortPadicMode); + MUNLOCK(ErrorMessageLock); + Terminate(-1); + return(0); + } + + fun3 = TermMalloc("MergeWithPadic"); + if ( AT.SortPadicMode == 3 ) padicret = AddPadics(BHEAD fun3,fun1,fun2); + else if ( AT.SortPadicMode == 1 ) + padicret = AddRatToPadic(BHEAD fun3,fun1,(UWORD *)coef2,size2); + else padicret = AddRatToPadic(BHEAD fun3,fun2,(UWORD *)coef1,size1); + if ( padicret != 0 ) { + TermFree(fun3,"MergeWithPadic"); + AT.SortPadicMode = 0; + return(0); + } + 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 : +*/ + +/* + #[ ConvertPadicLimbs32 : +*/ +int ConvertPadicLimbs32(int32_t *out, LONG capacity, const int32_t *in, + int sourcebits, int targetbits) +{ + LONG n, count, size, i; + uint32_t high; + if ( (sourcebits != 16 && sourcebits != 32) + || (targetbits != 16 && targetbits != 32) || in[0] != PADICLIMBS + || in[2] == INT32_MIN || in[2] == 0 ) return(-1); + n = in[2] < 0 ? -(LONG)in[2] : in[2]; + if ( in[1] < PADICLIMBHEAD+1 || n != (LONG)in[1]-PADICLIMBHEAD + || in[1] > (sourcebits == 16 ? INT16_MAX : INT32_MAX) ) return(-1); + if ( sourcebits > targetbits && n > (INT16_MAX-PADICLIMBHEAD+1)/2 ) + return(-2); + high = (uint32_t)in[PADICLIMBHEAD+n-1]; + if ( sourcebits == 16 ) high &= 0xffffU; + if ( high == 0 ) return(-1); + count = n; + if ( sourcebits < targetbits ) count = (n+1)/2; + else if ( sourcebits > targetbits ) count = 2*n-(high <= 0xffffU); + size = PADICLIMBHEAD+count; + if ( size > capacity || size > (targetbits == 16 ? INT16_MAX : INT32_MAX) ) + return(-2); + if ( out ) { + out[0] = PADICLIMBS; out[1] = (int32_t)size; + out[2] = (int32_t)(in[2] < 0 ? -count : count); + for ( i = 0; i < count; i++ ) { + uint32_t limb; + if ( sourcebits < targetbits ) { + limb = (uint32_t)in[PADICLIMBHEAD+2*i] & 0xffffU; + if ( 2*i+1 < n ) limb |= ((uint32_t)in[PADICLIMBHEAD+2*i+1] & 0xffffU) << 16; + } + else if ( sourcebits > targetbits ) + limb = ((uint32_t)in[PADICLIMBHEAD+i/2] >> (16*(i%2))) & 0xffffU; + else limb = (uint32_t)in[PADICLIMBHEAD+i]; + out[PADICLIMBHEAD+i] = (int32_t)limb; + } + } + return((int)size); +} +/* + #] ConvertPadicLimbs32 : +*/ diff --git a/sources/pre.c b/sources/pre.c index 94c36332e..f6c7ea483 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} @@ -7786,6 +7792,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(); @@ -7888,5 +7900,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 b0fc8f74c..b6e97d1da 100644 --- a/sources/proces.c +++ b/sources/proces.c @@ -4007,6 +4007,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/reshuf.c b/sources/reshuf.c index 1836f3108..dd44149bb 100644 --- a/sources/reshuf.c +++ b/sources/reshuf.c @@ -141,6 +141,9 @@ void FunLevel(PHEAD WORD *term) switch ( *t ) { case SYMBOL: case DOTPRODUCT: +#ifdef WITHFLINT + case PADICLIMBS: +#endif break; case VECTOR: t += 3; diff --git a/sources/sch.c b/sources/sch.c index 24476f0e3..2f7e78457 100644 --- a/sources/sch.c +++ b/sources/sch.c @@ -1412,6 +1412,17 @@ void WriteArgument(WORD *t) lowestlevel = 0; if ( *t > 0 ) { oldoutsidefun = AC.outsidefun; AC.outsidefun = 0; +#ifdef WITHFLINT + if ( t[ARGHEAD] == *t-ARGHEAD && TestPadicLimbTerm(t+ARGHEAD) ) { + LONG size = ARGHEAD+2*(LONG)ABS(t[ARGHEAD+PADICLIMBHEAD])+2; + WORD *expanded; + expanded = Malloc1((size_t)size*sizeof(WORD),"PrintPadicUnit"); + ExpandPadicArgument(expanded,t); + WriteExpression(expanded+ARGHEAD,size-ARGHEAD); + M_free(expanded,"PrintPadicUnit"); + } + else +#endif WriteExpression(t+ARGHEAD,(LONG)(*t-ARGHEAD)); AC.outsidefun = oldoutsidefun; goto CleanUp; @@ -1939,6 +1950,11 @@ int WriteInnerTerm(WORD *term, WORD first) WORD *t, *s, *s1, *s2, n, i, pow; #ifdef WITHFLOAT int FloatChars = 0; +#endif +#ifdef WITHFLINT + WORD *PrintedPadic = 0; +#endif +#ifdef WITHFLOAT GETIDENTITY #endif t = term; @@ -2016,6 +2032,28 @@ 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; + first = 1; + while ( ss < t ) { + if ( *ss == PADICFUN ) { + if ( PrintPadic(ss) != 0 ) { + TokenToLine(AO.padicspace); + PrintedPadic = ss; + first = 0; + break; + } + } + ss += ss[1]; + } + } #endif else first = 1; while ( s < t ) { @@ -2066,6 +2104,11 @@ int WriteInnerTerm(WORD *term, WORD first) if ( *s == FLOATFUN && AO.FloatPrec >= 0 && AT.aux_ != 0 ) { } else +#endif +#ifdef WITHFLINT + if ( s == PrintedPadic ) { + } + else #endif { if ( *s >= FUNCTION && AC.funpowers > 0 diff --git a/sources/sort.c b/sources/sort.c index 69de71e76..dcfcc65e6 100644 --- a/sources/sort.c +++ b/sources/sort.c @@ -1802,6 +1802,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; @@ -2419,6 +2422,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); @@ -2427,6 +2433,16 @@ WORD Compare1(PHEAD WORD *term1, WORD *term2, WORD level) t1 = term1 + 1; t2 = term2 + 1; while ( t1 < stopper1 && t2 < stopper2 ) { +#ifdef WITHFLINT + if ( level == 0 && AC.activePadic ) { + WORD padic1 = *t1 == PADICFUN && t1+t1[1] == stopper1 && TestPadic(t1); + WORD padic2 = *t2 == PADICFUN && t2+t2[1] == stopper2 && TestPadic(t2); + /* A terminal p-adic coefficient ends the monomial, just like a rational. */ + if ( padic1 && padic2 ) { AT.SortPadicMode = 3; return(0); } + if ( padic1 ) return(PREV(AR.SortType == SORTLOWFIRST ? 1 : -1)); + if ( padic2 ) return(PREV(AR.SortType == SORTLOWFIRST ? -1 : 1)); + } +#endif if ( *t1 != *t2 ) { if ( *t1 == HAAKJE ) return(PREV(-1)); if ( *t2 == HAAKJE ) return(PREV(1)); @@ -2458,6 +2474,18 @@ WORD Compare1(PHEAD WORD *term1, WORD *term2, WORD level) } if ( c1 <= (FUNCTION-1) || ( c1 >= FUNCTION && functions[c1-FUNCTION].spec > 0 ) ) { +#ifdef WITHFLINT + if ( c1 == PADICLIMBS ) { + /* Limbs span the full unsigned WORD range; never subtract them. */ + while ( s1 < t1 && s2 < t2 ) { + if ( *s1 != *s2 ) return(PREV(*s1 < *s2 ? 1 : -1)); + s1++; s2++; + } + if ( s1 < t1 ) return(PREV(-1)); + if ( s2 < t2 ) return(PREV(1)); + } + else +#endif if ( c1 == SYMBOL ) { if ( *s1 == FACTORSYMBOL && *s2 == FACTORSYMBOL && s1[-1] == 4 && s2[-1] == 4 @@ -2819,6 +2847,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 ) { @@ -3992,6 +4030,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 126b7d3b1..170023ac9 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/store.c b/sources/store.c index cfa21392d..91605e808 100644 --- a/sources/store.c +++ b/sources/store.c @@ -1990,6 +1990,9 @@ void DetVars(WORD *term, WORD par) term++; } } +#ifdef WITHFLINT + else if ( *term == PADICLIMBS ) term = t; +#endif else if ( *term == HAAKJE ) term = t; else { if ( *term > MAXBUILTINFUNCTION ) @@ -2511,6 +2514,9 @@ static LONG ctrap=0; /*!!! ctrap++; !!!*/ +#ifdef WITHFLINT + if ( *term == PADICLIMBS ) { term += term[1]; continue; } +#endif if ( *term == SYMBOL ) { t = term + term[1]; term += 2; @@ -3967,42 +3973,42 @@ static void ResizeCoeff32(UBYTE **bout, UBYTE *bend, UBYTE *top) } else { - /* 2 -> 4 */ + /* 2 -> 4. The 16-bit limbs are stored in signed 32-bit cells. */ int32_t len = (end - 1 - out) / 2; if ( len == 1 ) { - *out = *(uint16_t *)out; + *out = (uint16_t)*out; ++out; - *out = *(uint16_t *)out; + *out = (uint16_t)*out; ++out; ++out; } else { p = out; - *out = *(uint16_t *)out; + *out = (uint16_t)*out; in = out + 1; for ( i = 1; i < len; ++i ) { /* shift */ - *out = (uint32_t)(*(uint16_t *)out) - + ((uint32_t)(*(uint16_t *)in) << 16); + *out = (uint32_t)((uint16_t)*out) + + ((uint32_t)((uint16_t)*in) << 16); ++in; if ( ++i == len ) break; /* copy */ ++out; - *out = *(uint16_t *)in; + *out = (uint16_t)*in; ++in; } ++out; - *out = *(uint16_t *)in; + *out = (uint16_t)*in; ++in; for ( i = 1; i < len; ++i ) { /* shift */ - *out = (uint32_t)(*(uint16_t *)out) - + ((uint32_t)(*(uint16_t *)in) << 16); + *out = (uint32_t)((uint16_t)*out) + + ((uint32_t)((uint16_t)*in) << 16); ++in; if ( ++i == len ) break; /* copy */ ++out; - *out = *(uint16_t *)in; + *out = (uint16_t)*in; ++in; } ++out; @@ -4909,6 +4915,32 @@ ReadSaveTerm32(UBYTE *bin, UBYTE *binend, UBYTE **bout, UBYTE *boutend, UBYTE *t newsubtermp = out; ++out; ++r; +#ifdef WITHFLINT + if ( id == PADICLIMBS ) { + int32_t *object = newsubtermp-1, *converted; + LONG oldsize = *newsubtermp, delta, limit = AM.MaxTer/sizeof(WORD); + int newsize; + if ( oldsize < PADICLIMBHEAD+1 || oldsize > coeff-object ) return(NULL); + newsize = ConvertPadicLimbs32(NULL,AM.MaxTal+PADICLIMBHEAD,object, + 8*AO.SaveHeader.lenWORD,BITSINWORD); + if ( newsize < 0 ) return(NULL); + delta = newsize-oldsize; + if ( end-newtermsize+delta > limit + || end-newtermsize+delta > WORD_MAX_VALUE ) return(NULL); + if ( (UBYTE *)(end+delta) > top ) return(bin); + converted = Malloc1((size_t)newsize*sizeof(int32_t),"LoadPadicLimbs"); + ConvertPadicLimbs32(converted,newsize,object, + 8*AO.SaveHeader.lenWORD,BITSINWORD); + r += oldsize-2; + out = object+newsize; + /* The caller moves the enclosing tail after a recursive resize. */ + memmove(out,r,(size_t)(end-t)*sizeof(int32_t)); + memcpy(object,converted,(size_t)newsize*sizeof(int32_t)); + M_free(converted,"LoadPadicLimbs"); + coeff += delta; end += delta; t += delta; + } + else +#endif if ( id == SYMBOL ) { while ( out < t ) { ++out; ++r; /* symbol number */ @@ -5031,7 +5063,12 @@ ReadSaveTerm32(UBYTE *bin, UBYTE *binend, UBYTE **bout, UBYTE *boutend, UBYTE *t r += t - out; out = t; } - else if ( id == GAMMA || id == LEVICIVITA || (id >= FUNCTION && AO.tensorList[id]) ) { + else if ( id == GAMMA || id == LEVICIVITA + || (id >= FUNCTION +#ifdef WITHFLINT + && id != PADICFUN +#endif + && AO.tensorList[id]) ) { /* int32_t vectoroffset = -2 * *((int32_t *)AO.SaveHeader.wildoffset); */ void *dummy = (void *)AO.SaveHeader.wildoffset; /* to remove a warning about strict-aliasing rules in gcc */ int32_t vectoroffset = -2 * *((int32_t *)dummy); @@ -5132,6 +5169,9 @@ ReadSaveTerm32(UBYTE *bin, UBYTE *binend, UBYTE **bout, UBYTE *boutend, UBYTE *t int32_t **ppp = &out; /* to avoid a compiler warning */ /* recursion */ newbin = ReadSaveTerm32((UBYTE *)r, binend, (UBYTE **)ppp, (UBYTE *)rend, top, 1); +#ifdef WITHFLINT + if ( newbin == NULL ) return(NULL); +#endif r += lenbuf; if ( newbin == (UBYTE *)r ) { return ( (UBYTE *)in ); @@ -5282,6 +5322,9 @@ LONG ReadSaveExpression(UBYTE *buffer, UBYTE *top, LONG *size, LONG *outsize) outend = out+*size; if ( outend > top ) outend = top; p = ReadSaveTerm32(in, end, &out, outend, top, 0); +#ifdef WITHFLINT + if ( p == NULL ) return(-1); +#endif if ( p == in ) break; in = p; } while ( in <= end - lenW ); @@ -5311,6 +5354,9 @@ LONG ReadSaveExpression(UBYTE *buffer, UBYTE *top, LONG *size, LONG *outsize) /* then renumber */ do { p = ReadSaveTerm32(in, end, &out, buffer+half, buffer+half, 0); +#ifdef WITHFLINT + if ( p == NULL ) return(-1); +#endif if ( p == in ) break; in = p; } while ( in <= end - sizeof(WORD) ); diff --git a/sources/structs.h b/sources/structs.h index 707e3e0f1..a18ae1a2e 100644 --- a/sources/structs.h +++ b/sources/structs.h @@ -1783,6 +1783,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 @@ -2113,6 +2118,9 @@ struct T_const { void *mpf_tab2; void *aux_; void *auxr_; +#endif +#ifdef WITHFLINT + void *padic_aux_; #endif NORMDATA **NormData; LONG NormDataSize; @@ -2185,6 +2193,10 @@ struct T_const { WORD FloatPos; WORD SortFloatMode; #endif +#ifdef WITHFLINT + WORD PadicPos; + WORD SortPadicMode; +#endif }; /* #] T : @@ -2380,6 +2392,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. */ @@ -2407,6 +2424,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 9dc5760c0..d9fbdf51f 100644 --- a/sources/threads.c +++ b/sources/threads.c @@ -4115,6 +4115,13 @@ int SortBotMerge(PHEAD0) *term1 = fun1-term1; TermFree(fun3,"SortBotMasterMerge"); } +#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 ); @@ -4773,4 +4780,4 @@ extern void optimize_expression_given_Horner_threaded(void) { #] optimize_expression_given_Horner_threaded : */ -#endif \ No newline at end of file +#endif diff --git a/sources/tools.c b/sources/tools.c index 65d3df532..1c1554fc7 100644 --- a/sources/tools.c +++ b/sources/tools.c @@ -3837,6 +3837,11 @@ int TestTerm(WORD *term) tstop = endterm - coeffsize; while ( t < tstop ) { switch ( *t ) { +#ifdef WITHFLINT + case PADICLIMBS: + if ( t != term+1 || !TestPadicLimbTerm(term) ) goto subtermsize; + break; +#endif case SYMBOL: case DOTPRODUCT: case INDEX: diff --git a/sources/transform.c b/sources/transform.c index 2ae4c4771..645961b95 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;