On Sun, May 12, 2024 at 04:43:53PM +0200, Grégory Vanuxem wrote:
> Hello,
>
> Le sam. 11 mai 2024 à 21:25, Waldek Hebisch <de...@fricas.org> a écrit :
> >
> > On Sat, May 04, 2024 at 09:30:51AM +0200, Grégory Vanuxem wrote:
> > > As a matter of fact, use of FLINT in FriCAS:
> > >
> > > (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
> >
> > 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).
>
> 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.
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
--
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/fricas-devel/ZkD9tO-NIrXGSbeS%40fricas.org.
PA := PrimitiveArray(INT)
Up := UnivariatePolynomial('x, INT)
do_unpack(i : INT, k : INT, l0 : INT, n0 : INT, mask : INT, ra : PA) : Void ==
n0 > 5 =>
n1 := n0 quo 2
i1 := shift(i, -k*n1)
do_unpack(i - shift(i1, k*n1), k, l0, n1, mask, ra)
do_unpack(i1, k, l0 + n1, n0 - n1, mask, ra)
for j in 0..(n0 - 1) repeat
ra(l0 + j) := LOGAND(i, mask)$Lisp pretend INT
i := shift(i, -k)
do_corr(t : INT, m : INT, n : INT, ra : PA) : Void ==
carry : INT := 0
for j in 0..(n - 1) repeat
i1 := ra(j) + carry
if i1 > t then
i1 := i1 - m
carry := 1
else
carry := 0
ra(j) := i1
unpack_int(i : INT, k : INT) : PA ==
l := length(i)
n := (l quo k)::NNI + 1
print [k, n]
ra := new(n, 0)$PA
if i < 0 then
i := i + shift(1, n*k)
do_unpack(i, k, 0, n, mask(k), ra)
m := shift(1, k)
do_corr(shift(m, -1) - 1, m, n, ra)
ra
pa_to_pol(pa : PA) : Up ==
n := (#pa - 1)::NNI
res : Up := 0
for j in 0..n repeat
res := monomial(pa(j), j)$Up + res
res
int_to_pol(i : INT, k : INT) : Up ==
pa := unpack_int(i, k)
pa_to_pol(pa)
