restart: interface(version);
theNumber:=100;
theFactorial:='theNumber!'; 
`~~`=evalf(%);;
                                             158
                     ~ = 0.933262154439442 10

# prime number decomposition, only the maximal factor of 2 is needed
ifactors(theFactorial)[2]: 
H:=sort(%): 
H[1]: 

a:='`^`'(%[1],%[2]);
b:='theFactorial/a'; #ifactor(%);

'type(b, posint)': 'type(b, odd)': '%%'=%%, '%'=%, 'a*b=theFactorial';

                            a := ^(2, 97)

                               theFactorial
                          b := ------------
                                    a

   type(b, posint) = true, type(b, odd) = true, a b = theFactorial


# Ansatz: 1) x = 0 modulo a and 2) x+1 = 0 modulo b, i.e. x = -1 = b-1 modulo b 
# this is a linear system
SYS:='s[1]*a = (b-1) + s[2]*b';
SOL:='isolve(%, t)';

                    SYS := s[1] a = b s[2] + b - 1

              SOL := isolve(s[1] a = b s[2] + b - 1, t)


eval(SYS, SOL): 
X:=lhs(%);       # parametric solution over the integers 
#'type(X, linear(t))': '%'=%; # it is linear in t


# take the positive value nearest to zero
t0:=ceil( solve(0=X) );
x:=eval(X, t=t0);

                               t0 := 0

  x := 16337039244185717783761368470132707882798968240021251359061\
        3415173983497717109834226880980122669349706207993240279765\
        8966752830270331977392256259918212890624

# this specific value solution of SYS also satisfies conditions by construction
'modp(x, a)': '%'=%;
'modp(x+1, b)': '%'=%;
                            modp(x, a) = 0
                          modp(x + 1, b) = 0

# thus x*(x+1) is divisible by both a and b:
'modp(x*(x+1), theFactorial)': '%'=%;

                  modp(x (x + 1), theFactorial) = 0

# now just verify desired constraints

After looking at Kitonum's result (in the last piece of code here) it seemed
to have x+1 = 0 modulo primary components for odd primes. The systems 'sys'
sets up a condition for that, giving a parametric solution in t (just think
of primary decomposing Z/(n!), i.e. fundamental theorem for abelian groups).

For p=2 one needs a different additional condition (x=0 modulo 2^(...) ).

Solving gives a linear form X in t. For n! = 100! one recognizes Kintonum's x.

Then give it a try at t=0 (one can do for 0 < X). Testing shows: that is ok.

That's the sketch after all failing approaches, and not very scientific ...
but computational very fast (and it seems to be correct)

@Markiyan Hirnyk 

Ok, that is code. But do you also know the reason for the task (and why one wants < 4*99!) ?

Besides how Maple handles these for me those things are quite different, even if they (here) look the same.

A vector is an element of a vector space (or more generally <blabla>), here we are in finite dimension

After choosing some base one can write v in X in coordinates (of the (orientated) base).

A matrix 'is' a linear map f: X --> Y between vector spaces (dim < oo), after choosing (orientated) bases.

But generally it makes no sense to 'add' f and v.

The suggestions means to interprete v as matrix and thus as map (using implicitly the base and having Y = underlying field)

Would you mind to provide the link and algo to shorten things? And why this is of interest - or is just for recreation? Is x known explicitly?

So users say "please set 'by time' as default" and Maple's answer is "do it yourself".

Very nice.

May I suggest that you drop some motivating words on your engineering examples (and that the user may wish to start the animation for the curve)?

Also to drop some words on the vector product (especially the notation A &x B) would be good.

Sometimes I like to write 'A.B': '%'=%; or 'A.B = B.A'; is(%); for a pretty printing.

Short answer: sign and Heaviside is 'the same', except how they behave in x=0
and 1 point does not matter for integrals (or generics for higher dimensions).

I got it by starring for a long time. But think that the correct way would be
to understand/examine your *initial* problem (without putting it into the form
of solving using integrals).

You can have my my ugly sheet(s), of course (in case: please provide a valid
adress beyond gmail or other snoopers, else I upload to webspace only)

The structure now is: if c < 0 then Int( c * Q(a,b,w)/q(b), 3-dim ), where the
c depneds linear in k_ii and the other variables. Which avoids lots of checks
needed for two Heaviside's and one signum

It also tells, that integration w.r.t. b should be done last (in theory)

Theoretical one/Maple can integrate Heaviside by parts. Leading to Dirac, which
can *NOT* be done numerical directly by Maple - only after finding 'zero curves'
Which is like (or close) to Green's theorem (but without parametrization).

But for your general parameters S_i you will want to use a non-linear solver,
being feed with the above:

One can *not* solve for k_i in a linear way (know all the integral form above).
That would mean, that the solution would no longer depend on the S_i. And I 
guess you already would have recognized and said that. 


But I can not really solve the task, as already said. 

Find (appended as link) a simplification of your integrals for better handling
and evaluation, having only 1 Heaviside (choose Heaviside(0) at your desire).

The integrands B_i come from your sheet, if you avoid numerical evaluation but
use symbolic solutions there (so 'convert/rational' is not needed). 

It should work as long as your S_i are 'generic' - by that I mean that they do
not enforce situations where the task in Heaviside(0) has dimensional problems
(you would need generic anyway in hoping for a solution).

I have not tried to separate that into 'coefficients' in k_i and solve for it
(because I do not see, that it would be a linear problem in the k_i, because
the Heaviside makes it non-linear - else it now could be solved numerically).

http://axelvogt.de/temp/MP_triple_integral%203dim%20system%203_1.pdf
http://axelvogt.de/temp/MP_triple_integral%203dim%20system%203_1.mws

@Dima 

It becomes too complicated for me having those parameters S1, S2, and S3
even if they are Reals (you mean that, yes?).

I only can suggest to try it on your own along the ideas for specific values
(where I am not sure about the symmetries for arbitrary S_i).

The only thing I can contribute is:

Your (1.-.2178*sin(b)) has to be seen as 1+1/2*(-1+1/(Pi^(1/2)))*sin(b)

Your 0.13808e-3 and 0.69039e-4 might be see as 1/7242 and 1/14484
so one is 2 times the other, up to numerical rounding errors.

Beyond that (trivial) remarks I have to give up.

@Carl Love yes, I used alias(delta(x) = delta*x)

Yes. But that is the problem for the 'machine' as well. Using enough digits (!)
one may have a chance using Plouffe's inverter (or variations)

http://isc.carma.newcastle.edu.au/
http://oldweb.cecm.sfu.ca/projects/ISC/ISCmain.html

In general I would not expect a system to respond by symbolic suggestion,
if not being feed by *many* constellations. And Maple covers some.

IIRC Plouffe once provided his full sources.

@litun 

then why not simply take the whole product and divide out the very factor ?

kmax:=10^4;
CodeTools:-Usage({seq(13^k    mod 5, k= 1..kmax)}):
CodeTools:-Usage({seq(13 &^ k mod 5, k= 1..kmax)}):

  memory used=44.77MiB, alloc change=0 bytes, cpu time=718.00ms, real time=718.00ms
  memory used=0.53MiB, alloc change=0 bytes, cpu time=15.00ms, real time=14.00ms

kmax:=10^5;
CodeTools:-Usage({seq(13^k    mod 5, k= 1..kmax)}):
CodeTools:-Usage({seq(13 &^ k mod 5, k= 1..kmax)}):

  memory used=10.08GiB, alloc change=0 bytes, cpu time=3.22m, real time=3.22m
  memory used=5.34MiB, alloc change=0 bytes, cpu time=141.00ms, real time=139.00ms

The improvement is great.


Now be unfair and     use little Fermat:

kmax:=10^5;
x:= modp(13,5);
CodeTools:-Usage({seq( x^modp(k, 5-1)    mod 5, k= 1..kmax)}):
CodeTools:-Usage({seq( x &^ modp(k, 5-1) mod 5, k= 1..kmax)}):

  memory used=5.15MiB, alloc change=0 bytes, cpu time=109.00ms, real time=109.00ms
  memory used=6.87MiB, alloc change=0 bytes, cpu time=110.00ms, real time=109.00ms


kmax:=10^6;
x:= modp(13,5):
CodeTools:-Usage({seq(13 &^ k         mod 5, k= 1..kmax)}):
CodeTools:-Usage({seq(x^modp(k, 5-1)  mod 5, k= 1..kmax)}):

  memory used=53.42MiB, alloc change=3.82MiB, cpu time=1.54s, real time=1.45s
  memory used=51.50MiB, alloc change=0 bytes, cpu time=936.00ms, real time=905.00ms

For kmax:=10^6 I get 

  memory used=53.42MiB, alloc change=3.82MiB, cpu time=1.54s, real time=1.45s
  memory used=51.50MiB, alloc change=0 bytes, cpu time=936.00ms, real time=905.00ms

and for kmax:=10^7

  memory used=0.52GiB, alloc change=324.66MiB, cpu time=16.62s, real time=15.60s
  memory used=0.50GiB, alloc change=38.15MiB, cpu time=10.09s, real time=9.52s



At kmax:=10^8 my system gives up, "Error, (in CodeTools:-Usage) object too large"
so I am not sure what is spent on computing, generating the sequence or finally
handling the set {seq}


But is a good improvement, though incomplete in my opinion.

why "because" - where to read up more for that ?