Two minor comments:

• As written, your linear equations all describe lines passing through the origin, and so they define only unbounded regions. Instead of , do you mean ?
• Disregarding the degenerate cases, where one or more lines are parallel or identical, there are not (8-1)! = 5040 intersection points between pairs of lines, but only binomial(8,2) = 28.

At present I have coded the following:

• Generating some coefficients (the factor 5 is present to get some interesting plot):
  

  n := 8:
  gen := rand(10000) - 5000:
  a := [seq(gen(),i=1..n)];
  b := [seq(gen(),i=1..n)];
  c := [5*seq(gen(),i=1..n)];
  

• Generating the equations for the lines:
  

  eqs := [seq(a[i]*x + b[i]*y = c[i],i = 1..n)];
  

• Plotting the lines:
  

  plot(map(solve,eqs,y),x = -10..10);
  

  giving, for instance,
  

But I have no clear idea yet as to how to obtain the coordinates of some points inside those regions. Maybe the fact (unless I am fundamentally mistaken) that all bounded regions are convex polygons is helpful: most probably there exist some powerful procedure in such a case. Hopefully, some hardcore mathematician, which I am not, can help us here.

That looks absolutely great. Thank you so much

Please ignore my other posting - I had not expected further reply

Suppose L is a list of expressions A[i]*x + B[i]*y - C[i], and S is a subset of {$1..n}. Then applyop(`*`,S,L,-1) multiplies the expressions numbered by members of S by -1. Make each of the results r into an inequality r < 0, and then you can call LinearUnivariateSystem on the result. Thus:

> L := [seq](A[i]*x + B[i]*y - C[i], i=1..8);
  P:= map( t ->  
     SolveTools[Inequality][LinearUnivariateSystem](
       map(`<`, applyop(`*`,t,L,-1), 0), x),
     combinat[powerset]({$1..8}));

But to actually get a point in each region, I think linear programming (LPSolve in the Optimization package) might be better. If a region given by a set of inequalities A[i]*x + B[i]*y - C[i] < 0 is nonempty and bounded, you can find a point in it by maximizing t such that A[i]*x + B[i]*y - C[i] + t <= 0. If the region is unbounded, add a constraint t <= M for some large M (say 10^6; this might have to be adjusted). So I might try something like this:

> getpt:= proc(T)
    local R, t;
    R:= convert(map(s -> (s + t <= 0), T),set) 
      union {t <= 10^6};
    R:= Optimization[LPSolve](t, R, maximize = true);
    if R[1] <= 0 then NULL 
    else subs(R[2],[x,y])
    fi;
  end;
  
  map(t -> getpt(applyop(`*`,t,L,-1)), 
        combinat[powerset]({$1..8}));

If I might just ask for a little more help. I am struggling to understand what is going on in the first section of your code Robert.

I have written it out in a little fuller form

L := [seq](A[i]*x + B[i]*y - C[i], i=1..8):
> CPS:=combinat[powerset]({$1..8}):
> P:= map( t -> ( map(`<`, applyop(`*`,t,L,-1) , 0) ) , CPS);

But that last line is eluding me. In particular the t in the applyop

Thank for all the help so far.

Thank you Robert that is most helpful.

I have also followed up the issue of number of regions produced by lines on a plane a
and an article can be found on
PlatonicDivisionsofSpace

For finding a point inside each bounded region, please see my blog post Regions bounded by lines.