Program overcounts or undercount points.

Question

I am completely new to maple, my professor showed me a couple of things and i am struggling to write a useful program.

First of all, I want all the square distances possible from points on a hexagonal grid to the origin AND how many points satisfy each distance (I need to pick the distance with most points) I use this: ( is there a better way to do that? Perhaps more efficient. )

In the end, it prints out the higher populated distance and I can hand pick it. However, later, I want to find the points that give me that highly populated distance so I use this:

The problem is, it gives me fewer points than the first counter claims there to be, so the first thing is overcounting something or this code is undercounting. I am not sure what the error is. Any help would be appriciated.

EDIT: I am sorry the photos are so small, there is no clear way to change maple into latex. To my knowledge.

Answer 1

This question on programming in Maple would be more appropriate on stackoverflow.com (or www.mapleprimes.com the Maplesoft user forum).

What is the maximal count you get when p and q range from just -300 to 300?

restart;

# This procedure will run faster under evalhf.
F:=proc(counter::Array(datatype=float[8]), n::integer)
  local p,q,dd;
  for p from -n to n do
    for q from -n to n do
      if irem(p+q,2)=0 then
        dd:=(p^2+3*q^2)/4;
        counter[(dd)]:=counter[(dd)]+1;
      end if;
    end do;
  end do:
  NULL;
end proc:

Now let's run it at for a smaller number of points.

# run without evalhf at smaller N
N:=300:
C:=Array(0..(2*N+1)^2,datatype=float[8]):

F(C, N):
max(C); # maximal count in Array
                          72.

max[index](C); # position at which maximal count occurs
                         12103

C[%]; # check
                          72.

Now compare when run under faster evalhf mode. The results concur.

# run under faster evalhf mode and compare
N:=300:
C:=Array(0..(2*N+1)^2,datatype=float[8]):

evalhf( F(C, N) ):
max(C); # maximal count in Array
                          72.

max[index](C); # position at which maximal count occurs
                         12103

C[%]; # check
                          72.

Now run the larger example, which takes a long time if done outside of evalhf mode.

# run larger example under evalhf
N:=3000:
C:=Array(0..(2*N+1)^2,datatype=float[8]):
evalhf( F(C, N) ):

max(C); # maximal count in Array
                          192.

max[index](C); # position at which maximal count occurs
                        1983163

C[%]; # check
                          192.

C[375193];  # claimed by OP as most frequent value
                          144.

As for your other piece of code which constructs the list of 2-valued points using repeated concatenation of a list... don't build it that way because it's very inefficient. Here's another way, where each "point" contains just the p and q value. Change the [p,q] to you formula for the "points" if you want.

maxindex:=max[index](C);
L:=[seq(seq(`if`((p^2+3*q^2)/4=maxindex,[p,q],NULL),
            q=-N..N),p=-N..N)]:

If you look at the final list L (as I make it directly above) you should be able to notice that it has some nice symmetry. Eg, [-2815,-53] is a point and so are [-2815,53], [2815,-53], and [2815,53]. If this is a programming exercise then you could consider whether you can construct that list L in 1/4 of the time, by computing over just one quadrant and then generating the rest of the points using those. When you have produced the full list you could apply your formula by using map on L.

[edited] Revision to demonstrate applying the formula to point [p,q].

restart;
H:=proc(counter::Array(datatype=float[8]), n::integer)
  local p,q,dd;
  for p from 1 to n do
    for q from 1 to n do
      if irem(p+q,2)=0 then
        dd:=(p^2+3*q^2)/4;
        counter[(dd)]:=counter[(dd)]+1;
      end if;
    end do;
  end do:
  NULL;
end proc:

N:=3000:
C:=Array(0..(2*N+1)^2,datatype=float[8]):
evalhf( H(C, N) ):

maxindex:=max[index](C);
L:=[seq(seq(`if`((p^2+3*q^2)/4=maxindex,
              op([[p,q],[-p,q],[p,-q],[-p,-q]]),NULL),
        q=1..N),p=1..N)]:

nops(L);

func:=pt->[pt[1]/2,pt[2]*sqrt(3)/2]:

L[13], func(L[13]); # applying formula to one point in L

final:=map(func, L):

plots:-pointplot(final);