devaluate(Domain) at beginning of new domain in Spad

42 views
Skip to first unread message

Grégory Vanuxem

unread,
May 4, 2024, 12:44:22 AM5/4/24
to fricas...@googlegroups.com
Hello,

I have found 'devaluate' in some Spad files and it can be handy. Some
domains are parameterized and it is not possible, or at least I don't
know how to do it, to use ' is ' instead of ' has ' (Domain vs.
Category) for them:

if R is Foo then
new(n) == bar(n)
else
new...

it is not possible to use:
if Foo(p) then

Or like #1 in "SubDomain(Integer, #1 >= 0)" used for NonNegativeInteger?

I don't know how. So I use 'devaluate' and a 'has' after the 'is', in
fact '=' here, see the string below, to let me gather information
about the domain, size(), precision() etc. (the parameters) via the
exported operations of the parent category. Think of ZMOD in a
container as an example.

But, this is not my question, sometimes using 'devaluate' breaks
completely the code after, it's still compilable but not executable.
Any ideas why?
R is a Ring.
===== spad ===============
NRing :String := string CAR((devaluate(R)$Lisp))$Lisp

Rep := SExpression
pprint := true

getind(m) ==> concat(["getindex(", "refs,_"", jlId m, "_")"])
=====
And later, the above macro is incorrectly treated (jlId is defined in
a category, and call the JLREFID *method*):

>> System error:
There is no applicable method for the generic function
#<STANDARD-GENERIC-FUNCTION BOOT::JLREFID (1)>
when called with arguments
(|NemoIntegerMod|).
See also:
The ANSI Standard, Section 7.6.6


A work around is:
================================================
NRing := CAR((devaluate(R)$Lisp))$Lisp pretend String
--NRing :String := string CAR((devaluate(R)$Lisp))$Lisp

Rep := SExpression
pprint := true

getind(m) ==> concat(["getindex(", "refs,_"", jlId m, "_")"])
================================================

And all becomes right... I have already encountered this issue before
but never found where it came from, now, yes!
And even a workaround. Before, I think it was adding some code between
the 'devaluate' and the macro that did the job.

- Greg

Grégory Vanuxem

unread,
May 4, 2024, 12:56:08 AM5/4/24
to fricas...@googlegroups.com
To add to this, the jlId method is normally called with a class
reference as argument:
#<JLREF zzModRingElem 864811387508848656 {100304E3B3}>

And not (|NemoIntegerMod|) so the generated code seems totally wrong.

Qian Yun

unread,
May 4, 2024, 1:05:01 AM5/4/24
to fricas...@googlegroups.com


On 5/4/24 12:43, Grégory Vanuxem wrote:
> Hello,
>
> I have found 'devaluate' in some Spad files and it can be handy. Some
> domains are parameterized and it is not possible, or at least I don't
> know how to do it, to use ' is ' instead of ' has ' (Domain vs.
> Category) for them:
>
> if R is Foo then
> new(n) == bar(n)
> else
> new...
>
> it is not possible to use:
> if Foo(p) then
>
> Or like #1 in "SubDomain(Integer, #1 >= 0)" used for NonNegativeInteger?

I'm not sure I understand your problem 100%, are you looking for
something like PrimeField in ffdoms.spad?

if not prime?(p)$IntegerPrimesPackage(Integer) then


Also it's difficult to know your problem with 'devaluate' in your
code snippet. Is there a minimal reproducible example?

- Qian

Grégory Vanuxem

unread,
May 4, 2024, 1:36:25 AM5/4/24
to fricas...@googlegroups.com
The first part of the mail about the use of devaluate is because, here, when constructing a matrix I need a _very_ specialized matrix type, FLINT[1] will be used so I can not use the generic Matrix type. As an illustration I use this later in this domain:

"devaluate" use is to define function(s), here the matrix constructor, in a specialized Matrix ring, so later with the Ring name as a string I can do:

else if NRing = "NemoPrimeField" and R has FiniteFieldCategory then
  new(rows, cols, a) ==
    jlref(concat(["matrix_space(GF(", string(size()$R), "),", string(rows), ",",
      string(cols),")(fill(", getind(a), ",(",
        string(rows), ",", string(cols),")))"]))
else if NRing is "NemoIntegerMod" and R has Finite then
  new(rows, cols, a) ==
    jlref(concat(["matrix_space(residue_ring(ZZ,", string(size()$R), ")[1]),", string(rows), ",",
      string(cols),")(fill(", getind(a), ",(",
        string(rows), ",", string(cols),")))"]))
        else if NRing = "NemoRealBall" and R has JuliaArbitraryPrecision then
          new(rows, cols, a) ==
            jlref(concat(["matrix_space(ArbField", string(precision()$R), "),", string(rows), ",",
              string(cols),")(fill(", getind(a), ",(",
                string(rows), ",", string(cols),")))"]))

"has Finite" or "has FiniteFieldCategory" let me call size() . And yes that is not beautiful code, the use of concatenated strings for constructions.

Not very easy to explain, I admit.




--
You received this message because you are subscribed to the Google Groups "FriCAS - computer algebra system" group.
To unsubscribe from this group and stop receiving emails from it, send an email to fricas-devel...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/fricas-devel/81462d8e-17a9-4d2f-8589-e3b76c049f39%40gmail.com.

Grégory Vanuxem

unread,
May 4, 2024, 3:31:32 AM5/4/24
to fricas...@googlegroups.com
As a matter of fact, use of FLINT in FriCAS:

(16) -> p:=(2*x+2*x^5+13*x^9)::UP("x",INT)

             9      5
   (16)  13 x  + 2 x  + 2 x
                                        Type: UnivariatePolynomial(x,Integer)
                                                                  Time: 0 sec
(17) -> p:=p^5;

                                        Type: UnivariatePolynomial(x,Integer)
                                                                  Time: 0 sec
(18) -> degree(p^750)

   (18)  33750
                                                        Type: PositiveInteger
                                     Time: 28.17 (EV) + 3.17 (GC) = 31.34 sec
(19) -> degree(p^750)

   (19)  33750
                                                        Type: PositiveInteger
                         Time: 0.02 (IN) + 31.17 (EV) + 3.06 (GC) = 34.25 sec
(20) -> degree(p^750)

   (20)  33750
                                                        Type: PositiveInteger
                                     Time: 29.72 (EV) + 2.93 (GC) = 32.65 sec
(21) -> )cl prop p
(21) -> x:=x::NUP(NINT,"x")

   (21)  x
                                Type: NemoUnivariatePolynomial(NemoInteger,x)
                                                                  Time: 0 sec
(22) -> p:=2*x+2*x^5+13*x^9

   (22)  13*x^9 + 2*x^5 + 2*x
                                Type: NemoUnivariatePolynomial(NemoInteger,x)
                                                                  Time: 0 sec
(23) -> p:=p^5;

                                Type: NemoUnivariatePolynomial(NemoInteger,x)
                                                                  Time: 0 sec
(24) -> degree(p^750)

   (24)  33750
                                                        Type: PositiveInteger
                                                   Time: 0.09 (EV) = 0.09 sec
(25) -> degree(p^750)

   (25)  33750
                                                        Type: PositiveInteger
                                                   Time: 0.06 (EV) = 0.06 sec
(26) -> degree(p^750)

   (26)  33750
                                                        Type: PositiveInteger
                                                   Time: 0.07 (EV) = 0.07 sec

Of course, different aims for FriCAS and Nemo/FLINT in Julia, even if now Python and Julia allow a more object oriented utilisation.

If you have Julia in your path, and Nemo installed in it and you want to give it a try, a simple:
git clone -b jlfricas --depth=1 https://github.com/gvanuxem/fricas.git

will  fetch the source. With sbcl just a './configure --enable-julia' should do the trick.
DISCLAIMER: draft code.

- Greg



Qian Yun

unread,
May 4, 2024, 3:58:15 AM5/4/24
to fricas...@googlegroups.com


On 5/4/24 15:30, Grégory Vanuxem wrote:
> As a matter of fact, use of FLINT in FriCAS:
>

My experience with FLINT is that is uses dense representation
for univariate polynomial.

That representation makes it fast for multiplication (when polynomials
are dense).

But I think very spare polynomial will choke it, for example

(x^10000000+1)^n (for n is some small integer).


====

So for FriCAS, lacking a dense representation makes it
performs very bad at certain tasks (for example, multiplication).

The karatsuba method implemented in
UnivariatePolynomialMultiplicationPackage is not working
at all because of the spare representation we have.

- Qian

Waldek Hebisch

unread,
May 4, 2024, 9:40:51 AM5/4/24
to fricas...@googlegroups.com
On Sat, May 04, 2024 at 06:43:42AM +0200, Grégory Vanuxem wrote:
> Hello,
>
> I have found 'devaluate' in some Spad files and it can be handy.
<snip>
> But, this is not my question, sometimes using 'devaluate' breaks
> completely the code after, it's still compilable but not executable.
> Any ideas why?
> R is a Ring.
> ===== spad ===============
> NRing :String := string CAR((devaluate(R)$Lisp))$Lisp
>
> Rep := SExpression
> pprint := true
>
> getind(m) ==> concat(["getindex(", "refs,_"", jlId m, "_")"])
> =====
> And later, the above macro is incorrectly treated (jlId is defined in
> a category, and call the JLREFID *method*):
>
> >> System error:
> There is no applicable method for the generic function
> #<STANDARD-GENERIC-FUNCTION BOOT::JLREFID (1)>
> when called with arguments
> (|NemoIntegerMod|).
> See also:
> The ANSI Standard, Section 7.6.6
>
>
> A work around is:
> ================================================
> NRing := CAR((devaluate(R)$Lisp))$Lisp pretend String
> --NRing :String := string CAR((devaluate(R)$Lisp))$Lisp

Robust use is

op_of_PS : Symbol := CAR(devaluate(PS)$Lisp)$Lisp

Usage like in 'fmtjfricas.spad':

n: String := string CAR(devaluate(f)$Lisp)$Lisp

is error prone: essentially Spad compiler can call _any_ function
'string' which returns a String. Namely, '$Lisp' disables normal
Spad typechecking and Spad compiler belives you that type is
right one. Without further restictions (like 'string$Symbl')
Spad will assume that call is OK. But the Lisp expression
returns Lisp symbol, so return value of Symbol is fine.
'pretend String' is wrong, Lisp symbol is _not_ a string.
Lisp symbol can be treated as SExpression, but 'string' in
SExpression requires Lisp string so will fail on Lisp symbols.
One you messed types 'sbcl' is free to generate wrong code.


Waldek Hebisch

Grégory Vanuxem

unread,
May 4, 2024, 2:43:22 PM5/4/24
to fricas...@googlegroups.com
I didn't notice your email, I just saw it. Modified accordingly,
that's better! Thanks.
> --
> You received this message because you are subscribed to the Google Groups "FriCAS - computer algebra system" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to fricas-devel...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/fricas-devel/ZjY64A3C5CtsmZJD%40fricas.org.

Grégory Vanuxem

unread,
May 6, 2024, 1:35:44 AM5/6/24
to fricas...@googlegroups.com
Hello Qian,

Le sam. 4 mai 2024 à 09:58, Qian Yun <oldk...@gmail.com> a écrit :
>
>
>
> On 5/4/24 15:30, Grégory Vanuxem wrote:
> > As a matter of fact, use of FLINT in FriCAS:
> >
>
> My experience with FLINT is that is uses dense representation
> for univariate polynomial.
>
> That representation makes it fast for multiplication (when polynomials
> are dense).
>
> But I think very spare polynomial will choke it, for example
>
> (x^10000000+1)^n (for n is some small integer).

Wow, very specific here. In the real world, is it of use?
Cryptography? Theoretical mathematics?
I looked at SAGE and they also use FLINT for generic cases like the
one I mentioned.

I will likely take a look at PARI/GP, NTL and Singular, just to inform
me, I don't really know some libraries/softwares that handle
specifically well sparse representations of polynomials, For now I
will have a look at SAGE with the keyword :sparse when instantiating a
polynomial ring. Thanks for this example.

- Greg


>
> ====
>
> So for FriCAS, lacking a dense representation makes it
> performs very bad at certain tasks (for example, multiplication).
>
> The karatsuba method implemented in
> UnivariatePolynomialMultiplicationPackage is not working
> at all because of the spare representation we have.
>
> - Qian
>
> --
> You received this message because you are subscribed to the Google Groups "FriCAS - computer algebra system" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to fricas-devel...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/fricas-devel/cbee840e-8104-4f39-839e-87c91c962cac%40gmail.com.

Qian Yun

unread,
May 6, 2024, 5:14:34 AM5/6/24
to fricas...@googlegroups.com


On 5/6/24 13:35, Grégory Vanuxem wrote:
> Hello Qian,
>
> Le sam. 4 mai 2024 à 09:58, Qian Yun <oldk...@gmail.com> a écrit :
>>
>>
>>
>> On 5/4/24 15:30, Grégory Vanuxem wrote:
>>> As a matter of fact, use of FLINT in FriCAS:
>>>
>>
>> My experience with FLINT is that is uses dense representation
>> for univariate polynomial.
>>
>> That representation makes it fast for multiplication (when polynomials
>> are dense).
>>
>> But I think very spare polynomial will choke it, for example
>>
>> (x^10000000+1)^n (for n is some small integer).
>
> Wow, very specific here. In the real world, is it of use?
> Cryptography? Theoretical mathematics?

Well, is your example ((2*x+2*x^5+13*x^9)^5)^750 useful in real world?
;-)

I think it's better to implement both and let user to choose dense
or sparse representation, or make an abstract interface that the
algorithm itself determines whether to use dense or sparse
representation.

> I looked at SAGE and they also use FLINT for generic cases like the
> one I mentioned.
>
> I will likely take a look at PARI/GP, NTL and Singular, just to inform
> me, I don't really know some libraries/softwares that handle
> specifically well sparse representations of polynomials, For now I
> will have a look at SAGE with the keyword :sparse when instantiating a
> polynomial ring. Thanks for this example.
>
> - Greg

FLINT uses sparse representation for multivariate polynomial.

- Qian

Waldek Hebisch

unread,
May 11, 2024, 3:25:18 PM5/11/24
to fricas...@googlegroups.com
On Sat, May 04, 2024 at 09:30:51AM +0200, Grégory Vanuxem wrote:
> As a matter of fact, use of FLINT in FriCAS:
>
Hmm, I wonder what you are computing here. Could you give result
(time and resulting number) for

reduce(_+, map(length, coefficients(p^750)))

assuming that 'length' and 'coefficients' are available for Nemo
integeres/polynomials. Also, I wonder if you can monitor actual CPU time
for this computation. I mean, if computation is delegated to a separate
process then FriCAS will not know about time spent in another process
and just report interface time. And one can split such computation
into multiple cores, that reduces real time (if you have enough
cores).

Also, is this time for multiplication or for powering? For
example, what happens if you do:

p5_f750 := p^750;

pp := p5_f750*p5_f750;


--
Waldek Hebisch

Grégory Vanuxem

unread,
May 11, 2024, 7:11:39 PM5/11/24
to fricas...@googlegroups.com
I am standing up... No usual opening hours

Just for information, there is a "lack" in FriCAC, pow is  implemented by default in a factored form of multiplication. here, this is directly with FLINT.
 

--
                              Waldek Hebisch


--
You received this message because you are subscribed to the Google Groups "FriCAS - computer algebra system" group.
To unsubscribe from this group and stop receiving emails from it, send an email to fricas-devel...@googlegroups.com.

Grégory Vanuxem

unread,
May 11, 2024, 8:00:02 PM5/11/24
to fricas...@googlegroups.com
(5) -> p5_f750 := p^750;

                                Type: NemoUnivariatePolynomial(NemoInteger,x)
                                                   Time: 0.15 (EV) = 0.15 sec
(6) -> pp := p5_f750*p5_f750;

                                Type: NemoUnivariatePolynomial(NemoInteger,x)
                                                   Time: 3.09 (EV) = 3.09 sec
(7) -> pp := p5_f750*p5_f750;

                                Type: NemoUnivariatePolynomial(NemoInteger,x)
                                                   Time: 2.81 (EV) = 2.81 sec
(8) -> pp := p5_f750*p5_f750;

                                Type: NemoUnivariatePolynomial(NemoInteger,x)
                                                   Time: 3.09 (EV) = 3.09 sec


For degree function, the chosen function differs. I have absolutely not implemented polynomial arithmetic, a TODO thing. 

Le sam. 11 mai 2024 à 21:25, Waldek Hebisch <de...@fricas.org> a écrit :
--
You received this message because you are subscribed to the Google Groups "FriCAS - computer algebra system" group.
To unsubscribe from this group and stop receiving emails from it, send an email to fricas-devel...@googlegroups.com.

Grégory Vanuxem

unread,
May 12, 2024, 10:44:33 AM5/12/24
to fricas...@googlegroups.com
Hello,

Le sam. 11 mai 2024 à 21:25, Waldek Hebisch <de...@fricas.org> a écrit :
>
Looking more carefully. I have not implemented coefficients since
coefficients in Nemo/AbstractAlgebra returns an iterator. But if you
want a "list" as in LISP you can use collect on this iterator. The
problem is that polynomial of degree, say 9, will return all
coefficients;:

julia> p
13*x^9 + 2*x^5 + 2*x

julia> collect(coefs)
10-element Vector{ZZRingElem}:
0
2
0
0
0
2
0
0
0
13

So in principle, in FriCAS, the default 'coefficients' operation is used:

(8) -> coefficients(p)

(8) [13, 2, 2]
Type: List(NemoInteger)

But I have also not implemented length for NemoInteger since Nemo does
not have it unless doing low level stuff I guess (it's a GMP mpz_t as
far as I know), so for your question (not first powered):

(6) -> p

(6) 13*x^9 + 2*x^5 + 2*x
Type: NemoUnivariatePolynomial(NemoInteger,x)
Time: 0 sec
(7) -> pp

9 5
(7) 13 y + 2 y + 2 y
Type: UnivariatePolynomial(y,Integer)
Time: 0 sec
(8) -> reduce(_+,map(length, coefficients(p^750)::List(INT)))

(8) 3672690
Type: PositiveInteger
Time: 0.30 (IN) + 1.47 (EV) = 1.78 sec
(9) -> reduce(_+,map(length, coefficients(pp^750)::List(INT)))

(9) 3672690
Type: PositiveInteger
Time: 0.20 (EV) + 0.01 (GC) = 0.22 sec

What you're asking is too costly for my actual implementation. I
should not have used the word "arithmetic" this morning, my mistake,
sorry, what I meant is that I have not implemented all
_primitive_routines, all polynomial manipulation/information routines.

With p^5^750, and the different coercions, FriCAS silently stop and
the process returns to the bash shell :-(

And for your question about the number of threads used the Nemo
documentation says:
-----------------------------------------------------
Experimental threading support for flint

Enabling a threaded version of flint can be done by setting the
environment variable NEMO_THREADED=1. To set the actual number of
threads, use Nemo.flint_set_num_threads($numberofthreads).

If you need more information, please let me know.
-----------------------------------------------------
https://nemocas.github.io/Nemo.jl/latest/#Experimental-threading-support-for-flint

Using it:

-- Execute a julia command directly in FriCAS and displa its returned value
(11) -> )jud Nemo.flint_set_num_threads(1)
To use threaded flint, julia has to be started with NEMO_THREADED=1
(11) ->

And time given by FriCAS seems good, that's extremely fast.

- Greg

Waldek Hebisch

unread,
May 12, 2024, 1:34:47 PM5/12/24
to fricas...@googlegroups.com
Well, 'coefficients' or coercion from NemoInteger to Integer should
be quite cheap. 'length' on Integer is _very_ cheap. I asked
for this to make sure that 'p^750' is actually computed. 'degree'
of 'p^750' is simply '750*degree(p)' so can be computed without
computing 'p^750'. Also, Julia folks put floating point computations
in various places. Correctly computing lenghts of coefficients
from floating point result is possible, but unlikely.

Using known "fast" algoritm FriCAS should be able to compute
'p^(5*750)' on your machine in about 1.5s. The algorithm
is:
1) convert polynomial to an integer by evaluation in x = 2^bl
where bl is bound of size of coefficients of the power
2) compute integer power
3) convert integer back to polynomial

When done literally, it would take probably about 6-7 seconds
of which 4 sec is for computing power and the rest is for
convertion. In principle convertion back should be very fast,
but using only operations available at Lisp level one needs
to do shifting and masking which creates a lot of intermediate
integers which are not needed later.

One can make it faster by noticing that coefficients of your
polynomial make arithmetic progression with step 4. So one can
reduce amount of work by factor of 4.

Using attached routines computation is:

)set messages time on
pt := 13*x^2 + 2*x + 2
pt5 := pt^5
bl := 15751
i5 := eval(pt5, x = 2^bl);
i5_750 := i5^750;
pt5_750 := int_to_pol(i5_750, bl);
p5_750 := multiplyExponents(pt5_750, 4)*monomial(1, 5*750)$UP('x, INT);
reduce(_+, map(length, coefficients(p5_750)))

You should see that only steps that take significant time are
computation of 'i5_750' (main cost) and 'int_to_pol' call.

Similar approach works for multiplication, but is such case one also
needs faster routine for converting polynomials to integers.

> I
> should not have used the word "arithmetic" this morning, my mistake,
> sorry, what I meant is that I have not implemented all
> _primitive_routines, all polynomial manipulation/information routines.
>
> With p^5^750, and the different coercions, FriCAS silently stop and
> the process returns to the bash shell :-(

Well, '5^750' is _much_ larger than '5*750', so that is huge computation
expected to overflow memory of any possible computer (unless implementaion
is lazy and do not compute actual result).

--
Waldek Hebisch
polm.input
Reply all
Reply to author
Forward
0 new messages