Advertisement
mate2code

Square roots of Gray code * bit-reversal

Mar 18th, 2014
930
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
MatLab 1.92 KB | None | 0 0
  1. % Matlab calculation of the OEIS sequences A239303 and A239304
  2. % The class bin is found here: http://pastebin.com/K02wGsae
  3.  
  4. SeqCV = sym(zeros(1,9870)) ;
  5. CountCV = 1 ;
  6.  
  7. SeqPerms = [] ;
  8.  
  9. for m=1:140  % because triangular(140)=9870, and OEIS b-files can have up to 10000 entries
  10.  
  11.     % create Mat (binary square matrix of order m)
  12.     if mod(m,2) == 0  % if m is even
  13.         MatTRBL = logical(eye(m/2)) ;
  14.         MatTL = fliplr(MatTRBL) ;
  15.         MatBR = [ false(1,m/2) ; [ false(m/2-1,1) fliplr(logical(eye(m/2-1))) ] ] ;
  16.         Mat = [ MatTL MatTRBL ; MatTRBL MatBR ] ;
  17.     else  % if m is odd
  18.        
  19.         MatTL = fliplr(logical(eye(ceil(m/2)))) ;
  20.         MatBR = fliplr(logical(eye(floor(m/2)))) ;
  21.         MatTR = [ fliplr(MatBR) ; false(1,floor(m/2)) ] ;
  22.         MatBL = MatTR' ;
  23.         Mat = [ MatTL MatTR ; MatBL MatBR ] ;
  24.     end
  25.    
  26.     % compression matrix to compression vector (interpret columns of Mat as binary numbers)
  27.     for n=1:m
  28.         SeqCV(CountCV) = sym(bin( find(Mat(:,n))-1 )) ;
  29.         CountCV = CountCV+1 ;
  30.     end
  31.  
  32.     % interpret Mat as graph and the graph as permutation
  33.     Start = find(diag(Mat)) ;  % first element of the permutation is the position of the only 1 on the diagonal of Mat
  34.     Perm = [ Start Start ] ;  % dirty trick 1
  35.     while 1
  36.         Last = Perm(end-1) ;
  37.         This = Perm(end) ;
  38.         Next = setdiff(  find(Mat(This,:)) ,  Last  ) ;
  39.         Perm = [ Perm Next ] ;
  40.         if isempty(Next)
  41.             break
  42.         end
  43.     end
  44.     Perm = Perm(2:end) ;  % dirty trick 2
  45.     SeqPerms = [ SeqPerms Perm ] ;
  46.  
  47. end
  48.        
  49. StrCV = '';
  50. for m=1:9870
  51.     StrCV = [ StrCV num2str(m) ' ' char(SeqCV(m)) '\n' ] ;
  52. end
  53. fid = fopen('BfileCV.txt','wt');
  54. fprintf(fid,StrCV);
  55. fclose(fid);
  56.  
  57. StrPerms = '';
  58. for m=1:9870
  59.     StrPerms = [ StrPerms num2str(m) ' ' num2str(SeqPerms(m)) '\n' ] ;
  60. end
  61. fid = fopen('BfilePerms.txt','wt');
  62. fprintf(fid,StrPerms);
  63. fclose(fid);
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement