home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / TestMatrix / magic_r < prev    next >
Encoding:
Text File  |  1994-12-20  |  2.3 KB  |  137 lines
  1. #
  2. # Magic square
  3. # magic() returns a square matrix whose elements
  4. # are 1 through N^2
  5. # Magic squares have equal row and column sums.
  6. #
  7. # Found this algorithm in a Fortran routine
  8. # on an old VAX at the University? Seems to
  9. # work quite well. Don't know who the author is though?
  10. #
  11.  
  12. magic = function ( N )
  13. {
  14.   a = zeros (N, N);
  15.    
  16.   if (int(mod(N,4)) != 0)
  17.   {
  18.     if (int(mod(N,2)) == 0) { m = N/2; }
  19.     if (int(mod(N,2)) != 0) { m = N; }
  20.  
  21.     # Odd order or upper corner of even order
  22.   
  23.     for (j in 1:m)
  24.     {
  25.       for (i in 1:m)
  26.       {
  27.         a[i;j] = 0;
  28.       }
  29.     }
  30.   
  31.     i = 1;
  32.     j = int((m+1)/2);
  33.     mm = m*m;
  34.  
  35.     for (k in 1:mm)
  36.     {
  37.       a[i;j] = k;
  38.       i1 = i - 1;
  39.       j1 = j + 1;
  40.       if (i1 < 1) { i1 = m; }
  41.       if (j1 > m) { j1 = 1; }
  42.       if(int(a[i1;j1]) != 0)
  43.       {
  44.         i1 = i + 1;
  45.         j1 = j;
  46.       }
  47.       i = i1;
  48.       j = j1;
  49.     }
  50.     if (int(mod(N, 2)) != 0) { return a; }
  51.   
  52.     # Rest of even order
  53.   
  54.     t = m*m;
  55.     for (i in 1:m)
  56.     {
  57.       for (j in 1:m)
  58.       {
  59.         im = i + m;
  60.         jm = j + m;
  61.         a[i;jm] = a[i;j] + 2*t;
  62.         a[im;j] = a[i;j] + 3*t;
  63.         a[im;jm] = a[i;j] + t;
  64.       }
  65.     }
  66.   
  67.     m1 = int((m-1)/2);
  68.     if (m1 == 0) { return a; }
  69.     for (j in 1:m1)
  70.     {
  71.       swap (a, m, 1, j, m+1, j);
  72.     }
  73.     m1 = int((m + 1)/2);
  74.     m2 = m1 + m;
  75.  
  76.     swap (a, 1, m1, 1, m2, 1);
  77.     swap (a, 1, m1, m1, m2, m1);
  78.  
  79.     m1 = N + 1 - int((m - 3)/2);
  80.     if (m1 > N) { return a; }
  81.  
  82.     for (j in m1:N)
  83.     {
  84.       swap (a, m, 1, j, m+1, j);
  85.     }
  86.     return a;
  87.  
  88.   else
  89.  
  90.     k = 1;
  91.     for (i in 1:N)
  92.     {
  93.       for (j in 1:N)
  94.       {
  95.         a[i;j] = k;
  96.         if (int(mod(i,4)/2) == int(mod(j,4)/2))
  97.         {
  98.           a[i;j] = N*N + 1 - k;
  99.         }
  100.         k = k + 1;
  101.       }
  102.     }
  103.   }
  104.   
  105.   return a;
  106. };
  107.  
  108. #
  109. # Swap elements in a matrix in column-major fashion.
  110. # At present their is no error checking.
  111. #
  112.  
  113. # A:    the input matirx
  114. # N:    number of elements to swap
  115. # R1:    the starting row number
  116. # C1:    the starting column number
  117. # R2:    the starting row number
  118. # C2:    the stating column number
  119. #
  120.  
  121. swap = function (A, N, R1, C1, R2, C2)
  122. {
  123.   local (i, nr, nc, tmp, s1, s2)
  124.  
  125.   nr = A.nr; nc = A.nc;
  126.   s1 = (C1-1)*nr + R1;
  127.   s2 = (C2-1)*nr + R2;
  128.  
  129.   for (i in 0:N-1)
  130.   {
  131.     tmp = A[s1+i];
  132.     A[s1+i] = A[s2+i];
  133.     A[s2+i] = tmp;
  134.   }
  135.   return 1;
  136. };
  137.