Items tagged with integer integer Tagged Items Feed

Maple needs better user-level facilities for doing linear algebra over finite fields, particularly the integers mod n. For example there is no good way to solve a linear system Ax=B when B is a matrix. Obviously the LinearAlgebra:-Modular package is very good at what it does. Why can't there be some nice non-programmer routines which call it ? One alternative to using the mod operator is to have all the commands in the main LinearAlgebra package accept an optional last argument for the characteristic. For example: LinearAlgebra:-GaussianElimination(A, n); Then in the GaussianElimination command you could do something like:
I've released a Maple implementation of the F4 algorithm for computing Groebner bases. You can download it from the Maple Application Center here, or from my personal webpage here. The code requires Maple 10. It's faster than the Groebner[Basis] command for total degree orders, and it can be run in non-commutative algebras too. The coefficients are restricted to rational numbers or the integers mod p with p
I am trying to solve (numerically) the following system of differential equations using dsolve (Maple 9): DGLS:=seq(seq(diff(B[i,k](t),t)=add(add(W(j,l,i,k)*B[j,l](t)-W(i,k,j,l)*B[i,k](t),l=1..N),j=1..N),k=1..N),i=1..N); The problem is posed by the coefficients W(i,k,j,l). I want to make the W(i,j,k,l) dependent on the unknown functions B[i,k](t). For that purpose I wrote a Maple program W := proc(i::integer, k::integer, j::integer, l::integer) which calculates a "transition probability" Wtrans. If Wtrans does not depend on the B[j,l] (e.g. Wtrans:=1/N^2;), everything works nice. However, if I try to use the B[i,k] somehow (e.g. if NKrit_jl > B[j,l] then Wtrans := 0.5 fi;), then error messages from the procedure W result like "Error, invalid terms in product".
Suppose you want to sort a list L ( of numbers ) and also determine order of elements as a list of indices corresponding to elements of original list L, that is, you want such a integer list "I" that "seq(L[I[j],j=1..nops(L))" to be equivalent to sorted list. This functionality is present in MATLAB in "[B,IX] = sort(...)" syntax and i come up with this problem while trying to convert a MATLAB function (GaussQuadratureWeights) to maple. The procedure described resolves problem using a few MAPLE commands including MAPLE's built-in sort function.
What is an efficient, robust, way to extract an integer coefficient from a single term? My first thought was using lcoeff, however, it doesn't work if the term contains constants (say Pi) or floats. Currently I'm using patmatch,
icoeff := proc(t)
local k,x,kx;
  return `if`(patmatch(t, k::'nonunit'(integer)*x::anything,kx)
end proc:
map(icoeff, [0, 1, -3.0, -3, -0.,  3.0*Pi, 4*Pi, -12/5*I]);
                           [0, 1, 1, -3, 1, 1, 4, 1]
Consider the task of sorting a list of complex floating-point numbers by magnitude. First Attempt The usual method to do this in Maple is with the sort procedure. By passing a boolean-valued function that computes then compares the magnitudes of two complex numbers, we can sort the list. The following procedure shows how this is accomplished.
sort1 := proc(L)
    return sort(L, proc(z1,z2) abs(z1) <= abs(z2) end proc);
end proc:
A disadvantage of this approach is that the absolute-value procedure is called twice every time a pair of numbers is compared. For a long list, the time spent in the absolute value routine dominates the computation time.
i have a list in the form of [[x,y,z]] where x,y,z are integers, i want to find all the elements of the form [0,y,z]. how do i do that?
Similarly to Binary Arithmetic, Octal Arithmetic can be done using following module,
Silvexil announced in his blog a Binary Arithmetic package. Here is my version of it.
Here is the simplest program for obtaining the set of prime numbers less than or equal to n,
It is not Eratosthenes Sieve though. The program ES below implements the Eratosthenes Sieve,

In this post, daniel asks about how to compute the list of primes less or equal to some integer n. This is often called the Sieve of Eratosthenes problem, but the term "Sieve of Eratosthenes" actually refers to a specific algorithm for solving this problem, not to the problem itself.

This problem interested me, so I tried a few different ways of performing the task in Maple. As with my previous blog post, this exercise should be seen as an attempt to explore some of the issues faced when trying to make Maple code run faster, not as an attempt to find the fastest all-time Maple implementation of an algorithm to solve this problem.

Here are seven different implementations:

Implementation 1

sp1 := n->select(isprime,{$1..n}):

This first was suggested by alec in response to daniel's post. It's by far the simplest code, and somewhat surprisingly turns out to be (almost) the fastest of the seven.

I had thought it wouldn't be, because the subexpression {$1..n}builds the expression sequence of all integers from 1 to n, and the code spends time checking the primality of all kinds of numbers which are obviously composite.

> time( sp1( 2^16 ) );

Implementation 2


This second attempt avoids the principal problem of sp1, which is the construction of that expression sequence mentioned previously. However, the evaluation of all those `if` conditionals is also expensive, so this approach is essentially the same as sp1in runtime.

> time( sp2( 2^16 ) );

Implementations 3 and 4

sp3 := proc(N)
    local S, n;
    S := {};
    n := 2;
    while n < N do
        S := S union {n};
        n := nextprime(n);
    end do;
end proc:
sp4 := proc(N)
    local S, n;
    S := {};
    n := prevprime(N);
    while n>2 do
        S := S union {n};
        n := prevprime(n);
    end do;
    S union {2}
end proc:

These two approaches are essentially the same: sp3 uses nextprime and ascends, while sp4 uses prevprimeand descends.

The advantage to this approach is in avoiding all the unnecessary primality tests done by sp1 and sp2. Unfortunately, this advantage is offset by the fact that they rely on augmenting the set Sincrementally, which causes them to be slower than either of the previous two.

> time( sp3( 2^16 ) );

> time( sp4( 2^16 ) );

Implementation 5

There's really only one basic Maple command dealing with primes that we haven't yet used: ithprime. However, to use this effectively we need to know how far to go: i.e. what is the maximal m such that pm ≤ n?

Well, that's equivalent to asking how many prime numbers there are between 1 and n, or equivalently, what's the value of π(n), where π(n) is the prime counting function. This is implemented in Maple as numtheory[pi], so we'll use that in our code.

sp5 := N->{seq(ithprime(i), i=1..numtheory:-pi(N))}:

Happily, this turns out to be much faster than anything we've seen yet:

> time( sp5( 2^16 ) );

Implementation 6

Just for fun, I wondered whether speed might be improved by approximating π(n) with the Prime Number Theorem using Li, the logarithmic integral. (There are lots of other, better, approximations that we could use, but I'm not going to bother with those here.)

The approximation isn't perfect: as the help page for Li says, π(1000)=168, but Li(1000) ≅ 178. So we'll use ithprimeto find the first Li(N) primes (rounded down), remove all primes > N, and add in any primes ≤ N that might be missing.

sp6 := proc(N)
    local M, S, n;
    # Approximate number of primes with Prime Number Theorem
    M := trunc(evalf(Li(N)));
    S := {seq(ithprime(i),i=1..M)};
    # Do some cleanup: since the theorem provides only an
    # approximation, we might have gone too far or not far enough.
    # Remove primes that are too big; add in any that are missing
    S := select( type, S, integer[2..N] );
    n := ithprime(M);
    while n < N do
        S := S union {n};
        n := nextprime(n);
    end do;
end proc:

This is faster than the first four attempts, but is not faster than sp5. It's unlikely to be, since numtheory[pi] is fairly fast. However, its speed could be improved with a closer approximation to π(n).

> time( sp6( 2^16 ) );

Implementation 7

This last attempt, to see what we get, is to cast away all of Maple's built-in tools, and actually implement a real Sieve of Eratosthenes. We dynamically build up a set of primes, then check each successive number for divisibility by each member of this set, making sure to use short-circuit evaluation so we don't waste time doing divisions once we know something's composite.

sp7 := proc(N)
    local S, n;
    S := {};
    for n from 2 to N do
        # if a prime p in S divides n, it's not prime
        if not ormap( p->irem(n, p)=0, S ) then
            S := S union {n};
        end if;
    end do;
end proc:

This approach is, as we might expect, a lot slower than anything we've tried thus far, because it redoes a lot of things that Maple already does quite quickly. For comparison, I've done it for inputs of 2^12 and 2^16so you can see the blow-up.

> time( sp7( 2^12 ) );
> time( sp7( 2^16 ) );

So, the conclusion is that sp5 is the fastest of the bunch. I'm not sure how far its runtime generalizes; it may be as fast as it is largely because it uses a precomputed list of primes. However, even for larger inputs I suspect it is probably still faster than any of the other approaches above, simply because it avoids the creation of large intermediate data structures or needless checks that most of the others perform.

Recently I was asked in a private email about the fastest way of calculating of binomial coefficients mod 2 in Maple. It shouldn't be a problem for anybody reading my assembly dll creation manual. Anyway, here is the assembly code,
With Maple being able to access dll functions, it is easy to produce a function giving the epoch from the date and time. I give an example here for doing that with Open Watcom compiler included in Maple's distribution for Windows.
It is quite frustrating how slow map or zip acts over rtables (examples below). I find it quite useful to write a separate procedure and use the new compiler abilities in Maple 10.

stack limit...

August 05 2005 jakubi 1299 Maple

Chi^2 calculations above some "size" or "complexity", using Maple 9.5 and Global Optimization Toolbox (GOT), may produce after some time of calculation error messages like:

"Execution stopped: stack limit reached.
The kernel has been shut down. Further computation cannot be performed."

Seeking workarounds, I have looked for information at ?kernelopts for kernelopts(stacklimit), but it was not very useful:

"Limits may be raised or lowered. Maple limits may not be raised above any system defined hard limits. "

First 19 20 21 22 Page 21 of 22