function [Y,M] = A337358(N) % MATLAB function for OEIS sequence A337358. % Given an upper limit N, returns 3 column matrix Y: % Y: 1st column is n values, % 2nd column is A335775 (real part), % 3rd column is A337358 (imaginary part). % also returns M, a matrix representation of the sequence. % Dimensions of M are the max row value X max column value. z=1; V=[]; for a = 2:N for b = 2:a if mod((sqrt(a^2+b^2)),1)==0 V=[V;a b sqrt(a^2+b^2)]; end end end B=sortrows(V,3); Y=[0 0 0]; M=zeros(1,1); r=1; c=1; for n=1:N P=[]; while B(z,3)==n P=[P;B(z,1:2)]; %possible moves z=z+1; end R=[]; %possible permutations for q=1:size(P,1) R=[R; r-P(q,1) c-P(q,2); r-P(q,1) c+P(q,2); r+P(q,1) c-P(q,2); r+P(q,1) c+P(q,2); r-P(q,2) c-P(q,1); r-P(q,2) c+P(q,1); r+P(q,2) c-P(q,1); r+P(q,2) c+P(q,1)]; end R=[R; r-n c; r+n c; r c-n; r c+n]; S=[]; for b=1:size(R,1) if all(R(b,:)>0) S=[S; R(b,:) norm(R(b,:))]; end end U=sortrows(S,[3 1]); M(U(1,1),U(1,2)) = n; r=U(1,1); c=U(1,2); Y=[Y; n c-1 r-1]; end end