The routine LinearAlgebra:-LA_Main:-LUDecomposition appears to select the first nonzero entry it can find (in the current leading column).

The title is the sound of my head hitting the desk.

The routine LinearAlgebra:-LA_Main:-LUDecomposition appears to select the first nonzero entry it can find (in the current leading column).

The title is the sound of my head hitting the desk.

Many "high level" commands: solve, simplify, normal, gcd, etc. will convert floats to rational numbers. It's not something I agree with, but sometimes there is no alternative.

This nifty construction does a brute force search, which is exponential time in the length of the modulus. That is not a very good way to do arithmetic. Imagine subtracting two numbers by searching for the possible result.

This nifty construction does a brute force search, which is exponential time in the length of the modulus. That is not a very good way to do arithmetic. Imagine subtracting two numbers by searching for the possible result.

Implement the extended Euclidean algorithm. Here is a version in C:

/* inverse of A mod p */
long modinv(long A, long p)
{
        long a, b, q, t, x, y;
        a = p;
        b = A;
        x = 1;
        y = 0;
        while (b != 0) {
                t = b;
                q = a/t;
                b = a - q*t;
                a = t;
                t = x;
                x = y - q*t;
                y = t;
        }
        return (y < 0) ? y+p : y;
}

Note that long is typically a 32 bit signed integer, division (q = a/t) is truncating, and the last line returns y if y is greater than or equal to zero or y+p otherwise. The result is the inverse of A mod p. I don't know anything about Excel, so you're on your own from here.

Implement the extended Euclidean algorithm. Here is a version in C:

/* inverse of A mod p */
long modinv(long A, long p)
{
        long a, b, q, t, x, y;
        a = p;
        b = A;
        x = 1;
        y = 0;
        while (b != 0) {
                t = b;
                q = a/t;
                b = a - q*t;
                a = t;
                t = x;
                x = y - q*t;
                y = t;
        }
        return (y < 0) ? y+p : y;
}

Note that long is typically a 32 bit signed integer, division (q = a/t) is truncating, and the last line returns y if y is greater than or equal to zero or y+p otherwise. The result is the inverse of A mod p. I don't know anything about Excel, so you're on your own from here.

It's hard to help you without giving you the answer. You need to write a program. In Maple that would look like:

F := proc(n)

# program goes here

end proc;

The hash # is a Maple comment, so the line "program goes here" is not part of the program. Inside you will need some logic. That would look like:

if n=1 then
  # do something
elif n=2 then
  # do something else
else
  # if it's not the other cases do this
end if;

Programs can return values using the return statement:

return 1;

will return the value 1.

return x;

will return x, whatever that is. You might want to make x an input to your program, for example:

F := proc(n,x)

instead of the version above. Programs can call themselves recursively. The usual way to do this is with procname, which refers to the program in which it appears.

procname(2,x)

will call whatever procedure it appears in, with the arguments 2 and x. You should be able to write your program now.

I assume you mean floating point. As far as I know, none of the options in LinearSolve are likely to succeed, which is a very unfortunate and (I think) ridiculous situation. You can try method='SparseIterative' or method='SparseDirect', and the help page also recommends method='hybrid'. If those fail I would email your linear system to Maplesoft and complain that their expensive software can not solve your real world problem, what are you paying for, etc. You can also try making equations out of your system and calling fsolve. If you post the system online I'm sure some of us will take a crack at it. I for one am interested.

Maple's identify function uses the PSLQ algorithm to find an integer relation among a small set of basis numbers (e, pi, sqrt(2), etc) that is close to the given number. It is limited by the number of constants in the basis, which is actually quite small. Plouffe's inverter is a massive online database. I assume he has is own server side lookup code and the Maple code sends a request, but he also says you can download the whole thing and run it locally, which would explain the 2MB of Maple code. As for the relative safety of system calls and whatnot, the fact that source is available is more of a reassurance than what you get from most programs. Nobody is really immune to the trust issue, see http://cm.bell-labs.com/who/ken/trust.html

Maple input:

restart;
alias(epsilon=ep);
eqn := ep*x^4 - 3*x^2 -2*x + 1;
eqn1 := simplify(subs(x=ep^(n)*y,eqn)/epsilon^(1+4*n));
asym1 := sum(y[n]*ep^(n), n=0..3);
meth1 := subs(n=-1/2,y=asym1,eqn1);
step21 := collect(expand(meth1,ep=0),ep);
zeroth1 := subs(ep=0,step21);
soln01 := [solve(zeroth1=0,y[0])];
coeff(step21,epsilon,2);

The problem is due to fractional powers in coeff. Someone else may have a better suggestion, but you can try:

select(has, step21, epsilon^2);

What's a "monomial"? Do they mean "a single term in polynomial"? Basically yes, a monomial is a product of variables to non-negative powers, like x^2*y^3*z. Why this requirement for "monomials in a to divide monomials in f"? ... ANY term can be made to divide another just by putting one over the other with a "/" in between. No, divide means exact division with a quotient. For example, x^2*y divides x^2*y^3 with quotient y^2, whereas x^3 does not divide x^2*y. There is no monomial alpha such that x^2*y = alpha*x^3. You might propose alpha = y/x but that is not a monomial. The reason for all of these rules is to avoid infinite recursion. By rewriting terms only when exact division is possible, you are guaranteed to "bottom-out". Now in your case x is independent of a*b, so that is not really a problem. You might try applying algsubs to the numerator and denominator separately, or you could also try "simplification with side relations":

expr := c/(a*b) + a*b*d;
simplify(expr, {a*b=x});

This will operate on the numerator and denominator separately. That is, it will normalize the fraction to (c+a^2*b^2*d)/a/b (see the normal command) and then divide (with remainder) by a*b-x. This rewrites all multiples of a*b, replacing a*b by x, repeatedly, until no more multiples of a*b are present. Maple 11 also contains some more powerful algorithms for simplifying fractions with side relations, which are useful in more complicated situations.

What's a "monomial"? Do they mean "a single term in polynomial"? Basically yes, a monomial is a product of variables to non-negative powers, like x^2*y^3*z. Why this requirement for "monomials in a to divide monomials in f"? ... ANY term can be made to divide another just by putting one over the other with a "/" in between. No, divide means exact division with a quotient. For example, x^2*y divides x^2*y^3 with quotient y^2, whereas x^3 does not divide x^2*y. There is no monomial alpha such that x^2*y = alpha*x^3. You might propose alpha = y/x but that is not a monomial. The reason for all of these rules is to avoid infinite recursion. By rewriting terms only when exact division is possible, you are guaranteed to "bottom-out". Now in your case x is independent of a*b, so that is not really a problem. You might try applying algsubs to the numerator and denominator separately, or you could also try "simplification with side relations":

expr := c/(a*b) + a*b*d;
simplify(expr, {a*b=x});

This will operate on the numerator and denominator separately. That is, it will normalize the fraction to (c+a^2*b^2*d)/a/b (see the normal command) and then divide (with remainder) by a*b-x. This rewrites all multiples of a*b, replacing a*b by x, repeatedly, until no more multiples of a*b are present. Maple 11 also contains some more powerful algorithms for simplifying fractions with side relations, which are useful in more complicated situations.

I can't believe I didn't notice that - thanks. One cheap way to fix it is to use 2*B+1, ie: 2001.

...it doesn't work with the mix of positive and negative integers implied by the assignment to r... I think it works fine. What problem do you see ?