@KIRAN SAJJAN 

Before going further, read my first answer and correct your ic/bc for H(Theta) accordingly: as the corresponding ode is 

diff(H(Theta), Theta)-.5*h(Theta)+.5*H(Theta)

you need ONLY ONE condition, not two!

Is that clear?

Second point now: three points boundary problem is pure nonsense from a mathematical point of view.
If you want to find a function f(Theta) whose boundary conditions are  f(-Pi/4)=0 and f(Pi/4)=0 and which veryfies f(0)=1, f(0)=1 CANNOT BE A 3rd BOUNDARY CONDITION.
f(0)=1 is a constraint on f.
To find f(Theta) add for instance the boundary condition D(f)(-Pi/4)=something and solve an optimization problem such that for some value of something you get f(0)=1.

The procedire is described inthe attached file (Case M=Mvals[1] only, I suppose you will be capable to do the rest).

Download Ode_New_TWO_phase-2_mmcdara.mw

@KIRAN SAJJAN 

The error you get can be avoided ithis way (illustration for Ans[1] only)

Ans[1](-Pi/4)
Error, (in unknown) solution is only defined in the range HFloat(-0.785398163397447)..HFloat(0.0)

Ans[1](evalf(-Pi/4))
Error, (in unknown) solution is only defined in the range HFloat(-0.785398163397447)..HFloat(0.0)

LeftBound := evalf(-Pi/4) + 1e-8;
Ans[1](LeftBound)
                         -0.7853981535
[                                      
[Theta(-0.7853981535) = -0.7853981535, 
[                                      

  H(Theta)(-0.7853981535) = HFloat(0.9999999999999999), 

  f(Theta)(-0.7853981535) = HFloat(1.2939666992754303e-16), 

  /   d            \(-0.7853981535)                              
  |------- f(Theta)|                = HFloat(2.6146695267427543e\
  \ dTheta         /                                             

       /   d    /   d            \\                              
  -8), |------- |------- f(Theta)||(-0.7853981535) = HFloat(2.64\
       \ dTheta \ dTheta         //                              

  1761642312108), 

  g(Theta)(-0.7853981535) = HFloat(-3.57443270369908e-9), 

  /   d            \                                               
  |------- g(Theta)|(-0.7853981535) = HFloat(-0.3611469438350275), 
  \ dTheta         /                                               

  h(Theta)(-0.7853981535) = HFloat(1.0000000012643802), 

  /   d            \                                             
  |------- h(Theta)|(-0.7853981535) = HFloat(0.12774816275307546)
  \ dTheta         /                                             

  ]
  ]
  ]

The fact plot accepts the lft boundary -Pi/4 is likely due to the way the plotting algorithm determines the points, which seems to be accredited by the following commands

display(
   odeplot(Ans[1],[Theta,f(Theta)], Theta=-Pi/4..-Pi/4+1e-3, numpoints=1000, color=cols[1])
   , 'axes'= 'boxed'
   , labels=[Theta,'f']
  );
Warning, left boundary of range exceeds domain of the problem, it has been adjusted

# No warning here
display(
   odeplot(Ans[1],[Theta,f(Theta)], Theta=LeftBound..-Pi/4+1e-3, numpoints=1000, color=cols[1])
   , 'axes'= 'boxed'
   , labels=[Theta,'f']
  );

@dharr 

"Your C and C' can just be combined into a single combined C that I deal with - the counts for C={1}, C'={2,3,4} are the same as for C={1,2}, C'={3,4}, so I just count {1,2,3,4}."
Good catch!

The rest is just as good: briliant.

By the way my poor algorithm required 45s to get the result with L=5 and N=10^5 (cpu time = 2.7s for yours on my machine, which should be significantly reduced using @acer's ideas)

@FDS 

Why not converting in seconds or days?
It seems quite natural to me for Maple to wait for your choice before making the corresponding conversion, doesn't it?

Suppose you were handling yards, feet and inches: would you have asked why Maple doesn't automatically convert to inches, or feet?

@acer 

I'm going to do all these tests and see what is best for me.
Thanks again.

@acer 

Your first solution based on evalhf and @dharr's idea is indeed a great improvement of what I did.

Concerning the use of Compile the result seems biased because the generation of the compile code is very lengthy:

t0 := time[real]():
cpattern:=Compiler:-Compile(pattern, 'inmem'=false):
time[real]()-t0;

t0 := time[real]():
cpattern(N,nC,M8,C4,P4):
time[real]()-t0;
                             5.912
                             0.009

Nevertheless, as I said in my initial question and detailed in my last reply to @dharr, I need do to this kind of "counting" for a large number of values of C and P. So the compilation time should indeed become marginal at last.

@dharr 

Good idea! 

I thought about using some "embedded property" for instance:

• Let C and C' two sets of columns such that C is strictly included in C'.
   
• Let P a list of length |C| of  zeros and ones and P' a "completion" of P obtained by adding "in place" zeros or ones, for instance

  C  = [1, 3]
  P  = [0, 1]
  C' = [1, 2, 3, 5]
  P' = [0, 1, 1, 1]

• Let  M(P) is the submatrix of M whose columns C match the pattern P, then the submatrix M(P') of M whose columns C' match the pattern P', is a submatrix of M(P).
  With L=5 and M built as I said in my initial question, M(P) has about N/4 rows and thus the search of the columns C' which match the pattern P', ia about 4 times faster when done with M(P).

In some sense it is the Reverse_Idea.mw of yours


By the way, in case you would wonder where this problem comes from, here is little explanation.
The true problemis that one; I have  L detectors, each sending a signal ("1") when some event occurs (each detector is assigned to record a particular event).  
For any event E_i=1..L its probability of occurrence p_i is constant over the whole record time.
Then the record history of detector D_i can be modelled by a (stationary) Bernoulli random process with parameter p_i.
One question is: are the events E_i=1..L mutually independent or not?

Generally one uses a restricted notion of independence named pairwise independence.
Let X, Y, Z three random variables: saying that they are pairwise independent means that X and Y are independent, X and Z are independent and Y and Z are independent.

Mutual independence is a stronger constraint: it states that any subset C of the proper powerset Q of {X, Y, Z} is independent of any other subset C' of P if C and C' are disjoint.
Obviously mutual independece implies pairwise independence while the reciprocal is generally false.

The independence of C and C' means that Prob(C = P | C' = P') = Prob(C = P) where P and P' are any of all the possible outcomes of C and C' respectively.
As I deal with Bernoulli random variables, such outcomes are sequences of 0 and 1. For instance, if C refers to detectors 2 and 3, P is one of {[0, 0], [0, 1], [1, 0],[1, 1]}. The same way, if C' refers to, lets say K detectors (detectors 2 and 3 excluded), there exists 2^K possible outcomes for P'.
Checking for the independence of C and C' then means assessing if the 2^K+2 conditions Prob(C = P | C' = P') = Prob(C = P) can be reasonnably accepted.

The number of conditions to check in order to assess the mutual independence is equal to 232 for L=4 detectors, and to 1320 for L=5 detectors. This explains why one usually check pairwise independence (L*(L-1)/2 conditions to check).

@C_R 

Thanks for your reply.

@C_R 

I had started to write a reply before you published yours, but I posted it after yours.
Sorry for the inconvenience and feel free to delete mine.

@dharr 

It's a good way to react. I should follow your advice.

do you systematically refuse to upload your worksheet while using the big green arrow?

@JAMET 

No, I still don't understand but you don't help me while not giving substantial informations.
Has this something to do with what is persented here ? 

I still wonder what your last reply has to do with your original question

Nevertheless there is a piece of code which does three things:

1. Build a Pythagorean Triple (PT) by applying a random composition of the three matrices given in the above reference to the primitive PT < 3, 4, 5>.
    
2. Build the the ternary tree of depth N of all PTs constructed from  < 3, 4, 5> by applying at most (N-1) of these three same matrices.
    
3. Find what matrices were applied to get a known pythagorean triple.

PTs.mw

PS: as I use Maple 2015.2 the graphic is quite poor but more recent Maple's version will enable you to do more beautiful things.

@acer 

Thanks

@JAMET 

I still don't understand.: where those 6 triples come from?
What is the triple you are starting from and what is the relations between them?

Finally, what does this question have in common with your first one?

what do you want to achieve?
Q1: Do you want to find all Pythagorean Triples (PT) (a, b, c) such that a < A, b < B, c < C where (A, B, C) is a PT?

Q2: Is (A, B, C) a "reduced" PT in the sense that the common divisor of A, B, C is 1?

Here is procedure which finds (within one arbitrary permutation) all the PT (a, b, c) assuming the answer to question Q1 is YES and (A, B, C) is "reduced".

 

Download SPT.mw