@janhardo The number of base-10 digits of M_n = 2^n-1 (regardless of whether that's prime) can be quickly computed with logarithms without needing to actually compute M_n, like this

MersenneDigits:= (n::posint)-> ceil(evalf[2+ilog10(n)](n*log10(2))):
MersenneDigits(82589933);
                            24862048

Maple has no trouble computing M_n itself; however, I warn you not to display it onscreen!! Use a colon to suppress the output!

@lijr07 I think that this'll work in Maple 2019:

TypeTools:-RemoveType('Monomial'): #avoid silly warning
TypeTools:-AddType(
    'Monomial',
    proc(e)
    local Var:= 'Y'[integer $ 2], Pow:= {Var, Var^rational};
        e::Pow or 
        e::'specop'(
            {Pow, Not({constant, 'satisfies'(e-> hastype(e, Var))})},
            `*`
        )
    end proc
)
:

@janhardo 

The memory allocation required for my code is quite modest and is proportional to p (the exponent), i.e., it's O(p). Specifically, it's (3/4)*p bytes for the data (4 p-bit integers and 1 2*p-bit integer) plus a constant amount for Maple to load the procedures.

However, you did discover (and I can document) an extremely severe memory-leak bug in Maple's handling of very-large-integer arithmetic. Specifically, this bug occurs if and only if the prime p > 524221. That's far above the range that I tested. I found a surprisingly simple workaround: Just replace s^2 with `*`(s, s) (I think that that may require 1D input to not get automatically converted back to s^2). Here's the updated code, and I've made a few other very minor changes.

(*---
Compute N mod (2^n-1) via bitwise operations. N must be < 4^n, which'll
always be true if this is used in a repeated-squaring operation such as
in IsMersennePrimeExponent.
-----*)
BitMod:= (N::posint, n::posint)->
local B:= Bits:-Split(N, n);
    `if`(nops(B) < 2, `if`(ilog2(N+1) = n, 0, N), thisproc(add(B), n))
:
(*---
Check whether (2^p-1) is prime by Lucas-Lehmer algorithm.
-----*)
IsMersennePrimeExponent:= proc(p::And(prime, Not(2)))
local s:= 4;
    to p-2 do s:= BitMod(`*`(s,s), p) - 2 od;
    evalb(s=0)
end proc:
IsMersennePrimeExponent(2):= true
:

The next prime after 524221 is 524331. I just tested with this exponent using a total memory allocation of 146.1 Mb and a time of 32 minutes. The time required for the algorithm is essentially proportional to p^2 when the more-efficient large-integer multiplication kicks in. Specifically, it's O(p^2*ln(p)*ln(ln(p))). I suspect that the value of p that I just used is precisely where that more-efficient multiplication kicks in because that's the first prime value of p such that the GMP representation of 2^p-1 has more than 2^13 memory words. Extrapolating from that 32 minutes, I predict that your prime 1398269 would take 4 hours 12 minutes on my computer, an Intel Core I-7 @ 2.7 GHz. Extra processors have little effect on this (but can be used to process multiple primes simultaneously).

@lijr07 Thanks. I modified the Answer so that it now also handles inert exponents, as in x%^2.

@lijr07 What about sin(t+x)*x?

@lijr07 Do you want to consider t^(1/3)*Y[1,2] to be a monomial? How about t^sqrt(2)*x? And what about t^a*x where a is an unassigned symbol? And sin(t)*x? All these are possible if you let me know.

@nm Why don't you respond? This was an important Question. I know that you're on MaplePrimes, because you've asked two Questions since I posted this.

@janhardo You asked:

• L::list   = input parameter L in procedureand  must be a list as datatype

Correct.

• local S:= 0, i;  local variables S and  i used in procedure
  S is sum of list L and i is index varialble of do loop.

Correct.

• I is a index variable in the do loop : don't know if this i starts on 0 or 1 ?

The index of a for-do loop starts at 1 unless there is a from clause or an in clause. Example (from):
    for i from 2 to nops(L) do
In that case the index starts at 2. An in clause causes the elements of the container to be extracted without using explicit indices. For example,
    for x in L do
The value of x will start at the "first" element of L. I put "first" in quotes because for some containers other than lists it's not necessarily true that this element is L[1].

• nops(L) is the number of elements in the list : here two in L:=[1,2] 

Yes, that's true for lists and sets, but not for some other containers. The command numelems works for sets, lists, and some other containers.

• So i cannot start on 0 ( there is no 0 th element in List L ),

Correct; list indices always start at 1. Negative integers can be used to index a list from its end (i.e., backwards), but 0 is never a valid list index. 

• Don't know why i should use this : S:= S + L[i]  and why not S:=L[i]  ?

If you use S:= L[i], then there's no addition operator, so how could anything get added?

• od; is a shortcut for:  end do

Almost; od (no semicolon) is short for end do. Also fi is short for end if. There are no other similar shortcuts.

• S  stands here alone  only for ?

It is the return value of the procedure. That line could also be return S. If it's the last executed statement of the procedure, then the word return is not necessary.

@janhardo You wrote:

•  r+  is  summing up by 1 for a integer number (natural numbers)?

In the statement r+= x, the operator is +=; the operator is not +. This is an example of an updating-assignment operator. Here's how they work:

1. Let OP represent any of Maple's predefined binary infix (i.e., two operands with the operator between them) operators whose return value is of the same general type as its operands. So, OP is any of {+, -, *, /, ., ^, mod, ||, and, or, xor, implies, union, intersect, minus}. Let r be any assignable expression, i.e., anything that could be on the left side of :=. And let x be any object. Then
       r OP= x
   is equivalent to
       r:= r OP x.
2. The updating-assignment operators ++ and -- are unary (i.e., one-operand) operators akin to += and -= with the second operand being 1. So
       ++x
   is equivalent to 
       (x+= 1),
   and likewise for --.
3. Point (2.) shows the unary operator used in prefix form (i.e., appearing before the operand). These two operators can also be used in postfix form (i.e., appearing after the operand), in which case the unupdated value of the variable is returned, though the variable is still updated. For example, 
       x:= 1; print(x++, x);
   prints 1, 2; whereas
       x:= 1; print(++x, x);
   prints 2, 2.
4. There is one more updating assignment operator: ,=. It's a bit too complicated to explain all its nuances here. See the help page.
5. All of this is on the help page ?assignment.
6. These operations are not restricted to natural numbers. They will work for any operands for which the regular (non-updating) forms of the operators would work.

@vv Vote up. Here is my variation of your Wind. It returns the same results using different syntax.

Wind:= proc(L::list([realcons, realcons]), P::[realcons, realcons])
local  
    AT:= proc(p) option inline; arctan(p[2],p[1])/Pi end proc,
    S:= p-> [0, ()-> undefined, signum][2+signum(0, abs(p)-1, 0)](p)
; 
    add(S~((u-> u[..-2] - u[2..])(map(AT@`-`, [L[], L[1]], P))))
end proc
:

@ktoh You need to copy-and-paste the procedure from my Answer in addition to the code from Acer's Reply. Acer included my procedure in the Startup Region of his worksheet---a part that usually isn't visible on-screen---thus its presence isn't obvious.

@aksupriatna There is no special package needed. Please show your error message, and state what version you're using and whether you're using 2d or 1d input.

@sand15 The solution space for msolve(..., n) is the set of equivalence classes of remainders mod n, which is a finite set with n elements. When working with modular arithmetic at anything other than an elementary level, it goes without saying that objects that are written as ordinary integers are actually representatives of their equivalence classes. In other words, x means {x + k*n | k in Z}. This a totally standard notation in modular arithmetic. 

IMO, algsubs is junk, and it should be deprecated. I would only use it as a last resort. Indeed, I haven't used it in at least 16 years. At the very least, its help page should advise that there are several newer and better commands, such as simplify with side relations. 

As you said, you're using an approximation method, so of course it deviates from the exact solution. In particular, approximation methods for IVPs tend to deviate more the further you get from the initial point.

So, what kind of help are you expecting? There's nothing wrong with your code as far as I can see. The deviation is inherent to the method.