cosine

=-((((COS(PI()*(ROW())/(COLUMN()))^2-1))^2-epsilon)/ABS((((COS(PI()*(ROW())/(COLUMN()))^2-1))^2-epsilon))-1)/2 where epsilon is a small number. This gives the mod pattern also although not in a “clean” way.

Posted in Uncategorized | Comments Off on cosine

Convolution

These formulas produce table A051731 in the oeis.
http://oeis.org/A051731

Excel spreadsheet formula, American version:

=IF(COLUMN()=1, 1, IF(ROW()>=COLUMN(), SUM(INDIRECT(ADDRESS(ROW()-COLUMN()+1, COLUMN()-1, 4)&":"&ADDRESS(ROW()-1, COLUMN()-1, 4), 4))-SUM(INDIRECT(ADDRESS(ROW()-COLUMN()+1, COLUMN(), 4)&":"&ADDRESS(ROW()-1, COLUMN(), 4), 4)), 0))

Excel spreadsheet formula, European version:

=IF(COLUMN()=1; 1; IF(ROW()>=COLUMN(); SUM(INDIRECT(ADDRESS(ROW()-COLUMN()+1; COLUMN()-1; 4)&":"&ADDRESS(ROW()-1; COLUMN()-1; 4); 4))-SUM(INDIRECT(ADDRESS(ROW()-COLUMN()+1; COLUMN(); 4)&":"&ADDRESS(ROW()-1; COLUMN(); 4); 4)); 0))

Written differently:

T(n,1)=1, k>1: T(n,k) =  \sum\limits_{n=1}^{k-1} T(n-i,k-1) - \sum\limits_{n=1}^{k-1} T(n-i,k)

Mats Granvik (mats.granvik(AT)abo.fi)

Posted in Uncategorized | Comments Off on Convolution

Where the Möbius transform and matrix inversion agree.

//Program starts, written in Scilab - a Matlab clone
x=20; // size of matrix, don't change

T=zeros(x,x);

//a = fraction A001790/A046161

a=[1/1; 1/2; 3/8; 5/16; 35/128; 63/256;
231/1024; 429/2048; 6435/32768; 12155/65536;
46189/262144; 88179/524288; 676039/4194304;
1300075/8388608; 5014575/33554432;9694845/67108864;
300540195/2147483648; 583401555/4294967296;
2268783825/17179869184; 4418157975/34359738368]

for n=1:x;
T(n,1)=a(n); //set first column equal to fraction
end
T;

for n=2:x;
  for k=2:n;
    s = 0;
    for i=1:k-1;
      s=s+T(n-i,k-1)-T(n-i,k);//divisibility recurrence
    end
    T(n,k)=s;
  end
end
T;

U=inv(T); //Matrix inverse

//inverse Moebius transform of first column in matrix U
V=zeros(x,x);
for n=1:x;
  for k=1:n;
    if modulo(n,k)==0;
      V(n,k)=U(n/k,1);
      end
  end
end
b=sum(V,'c');
// a is the input and b is the output. We notice that a=b and
// therefore it could perhaps be said that the fraction
// A001790/A046161 is invariant under the divisibility recurrence,
// matrix inversion and the inverse Moebius transform.
C=[a,b]
//
// The divisibility recurrence applied to a sequence is the same
// thing as the Dirichlet convolution of another sequence.
Posted in Uncategorized | Comments Off on Where the Möbius transform and matrix inversion agree.

characteristic function of primes

-((A000005-2-1/2)/sqrt((A000005-2-1/2)*(A000005-2-1/2))-1)/2

Posted in Uncategorized | Comments Off on characteristic function of primes

Mertens function invariant

//Scilab (a Matlab clone)

x=10;  // size of matrix
T=zeros(x,x);

for n=1:x;
T(n,1)=1;
end
T;

for n=1:3;
  for k=1:n;
    T(n,k)=1;
  end
end
T;

// Random numbers here:
for n=4:x;
  T(n,1)=rand();
end
T;

for n=4:x;
  k=2;
  T(n,k)=T(n,k-1)-T(n-1,k);
end
T;

for n=4:x;
  k=3;
  T(n,k)=T(n,k-1)-(T(n-1,k)+T(n-2,k));
end
T;

for n=4:x;
  for k=4:n;
  s_one = 0;
    for i=1:k-2;
    s_one=s_one+T(n-i,k-1);
    end
  s_two = 0;
    for i=1:k-1;
    s_two=s_two+T(n-i,k);
    end
  T(n,k)=s_one-s_two;
  end
end
T

A = inv(T) // Mertens function in first columns of A
// Mats Granvik mats.granvik(AT)abo.fi
Posted in Uncategorized | Comments Off on Mertens function invariant

The inverse of a triangular matrix using determinants

Consider the lower triangular matrix A:
\left(\begin{array}{cccccc}a_{\mathrm{11}}&0&0&0&0&0\\ a_{\mathrm{21}}&a_{\mathrm{22}}&0&0&0&0\\ a_{\mathrm{31}}&a_{\mathrm{32}}&a_{\mathrm{33}}&0&0&0\\ a_{\mathrm{41}}&a_{\mathrm{42}}&a_{\mathrm{43}}&a_{\mathrm{44}}&0&0\\ a_{\mathrm{51}}&a_{\mathrm{52}}&a_{\mathrm{53}}&a_{\mathrm{54}}&a_{\mathrm{55}}&0\\ a_{\mathrm{61}}&a_{\mathrm{62}}&a_{\mathrm{63}}&a_{\mathrm{64}}&a_{\mathrm{65}}&a_{\mathrm{66}}\end{array}\right)

Modify as follows into matrix B:

\left(\begin{array}{cccccc}a_{\mathrm{11}}&0&0&0&0&1/(a_{\mathrm{11}}*a_{\mathrm{22}}*a_{\mathrm{33}}*a_{\mathrm{44}}*a_{\mathrm{55}}*a_{\mathrm{66}})\\ a_{\mathrm{21}}&a_{\mathrm{22}}&0&0&0&0\\ a_{\mathrm{31}}&a_{\mathrm{32}}&a_{\mathrm{33}}&0&0&0\\ a_{\mathrm{41}}&a_{\mathrm{42}}&a_{\mathrm{43}}&a_{\mathrm{44}}&0&0\\ a_{\mathrm{51}}&a_{\mathrm{52}}&a_{\mathrm{53}}&a_{\mathrm{54}}&a_{\mathrm{55}}&0\\ a_{\mathrm{61}}&a_{\mathrm{62}}&a_{\mathrm{63}}&a_{\mathrm{64}}&a_{\mathrm{65}}&0\end{array}\right)

Then calculate the determinant of matrix B and you will get the value at the sixth row and first column of the inverse of matrix A.

To calculate the value of the matrix inverse at the sixth row and second column calculate the determinant of the following matrix:

\left(\begin{array}{cccccc}a_{\mathrm{11}}&0&0&0&0&0\\ a_{\mathrm{21}}&a_{\mathrm{22}}&0&0&0&1/(a_{\mathrm{11}}*a_{\mathrm{22}}*a_{\mathrm{33}}*a_{\mathrm{44}}*a_{\mathrm{55}}*a_{\mathrm{66}})\\ a_{\mathrm{31}}&a_{\mathrm{32}}&a_{\mathrm{33}}&0&0&0\\ a_{\mathrm{41}}&a_{\mathrm{42}}&a_{\mathrm{43}}&a_{\mathrm{44}}&0&0\\ a_{\mathrm{51}}&a_{\mathrm{52}}&a_{\mathrm{53}}&a_{\mathrm{54}}&a_{\mathrm{55}}&0\\ a_{\mathrm{61}}&a_{\mathrm{62}}&a_{\mathrm{63}}&a_{\mathrm{64}}&a_{\mathrm{65}}&0\end{array}\right)

To calculate the value of the matrix inverse at the sixth row and third column calculate the determinant of the following matrix:

\left(\begin{array}{cccccc}a_{\mathrm{11}}&0&0&0&0&0\\ a_{\mathrm{21}}&a_{\mathrm{22}}&0&0&0&0\\ a_{\mathrm{31}}&a_{\mathrm{32}}&a_{\mathrm{33}}&0&0&1/(a_{\mathrm{11}}*a_{\mathrm{22}}*a_{\mathrm{33}}*a_{\mathrm{44}}*a_{\mathrm{55}}*a_{\mathrm{66}})\\ a_{\mathrm{41}}&a_{\mathrm{42}}&a_{\mathrm{43}}&a_{\mathrm{44}}&0&0\\ a_{\mathrm{51}}&a_{\mathrm{52}}&a_{\mathrm{53}}&a_{\mathrm{54}}&a_{\mathrm{55}}&0\\ a_{\mathrm{61}}&a_{\mathrm{62}}&a_{\mathrm{63}}&a_{\mathrm{64}}&a_{\mathrm{65}}&0\end{array}\right)

and so on.

See also: Inverse of a triangular matrix by matrix multiplication
mats.granvik(AT)abo.fi

Posted in Uncategorized | Comments Off on The inverse of a triangular matrix using determinants

The inverse of a triangular matrix by matrix multiplication

Here follows a technique for inverting triangular matrices that was discovered together with Gary W. Adamson. His blog can be found at http://qntmpkt.blogspot.com/

Consider the lower triangular matrix A:

\left(\begin{array}{cccccc}a_{\mathrm{11}}&0&0&0&0&0\\ a_{\mathrm{21}}&a_{\mathrm{22}}&0&0&0&0\\ a_{\mathrm{31}}&a_{\mathrm{32}}&a_{\mathrm{33}}&0&0&0\\ a_{\mathrm{41}}&a_{\mathrm{42}}&a_{\mathrm{43}}&a_{\mathrm{44}}&0&0\\ a_{\mathrm{51}}&a_{\mathrm{52}}&a_{\mathrm{53}}&a_{\mathrm{54}}&a_{\mathrm{55}}&0\\ a_{\mathrm{61}}&a_{\mathrm{62}}&a_{\mathrm{63}}&a_{\mathrm{64}}&a_{\mathrm{65}}&a_{\mathrm{66}}\end{array}\right)

Divide the COLUMNS with the diagonal elements in matrix A:

\left(\begin{array}{cccccc}a_{\mathrm{11}}/a_{\mathrm{11}}&0&0&0&0&0\\ a_{\mathrm{21}}/a_{\mathrm{11}}&a_{\mathrm{22}}/a_{\mathrm{22}}&0&0&0&0\\ a_{\mathrm{31}}/a_{\mathrm{11}}&a_{\mathrm{32}}/a_{\mathrm{22}}&a_{\mathrm{33}}/a_{\mathrm{33}}&0&0&0\\ a_{\mathrm{41}}/a_{\mathrm{11}}&a_{\mathrm{42}}/a_{\mathrm{22}}&a_{\mathrm{43}}/a_{\mathrm{33}}&a_{\mathrm{44}}/a_{\mathrm{44}}&0&0\\ a_{\mathrm{51}}/a_{\mathrm{11}}&a_{\mathrm{52}}/a_{\mathrm{22}}&a_{\mathrm{53}}/a_{\mathrm{33}}&a_{\mathrm{54}}/a_{\mathrm{44}}&a_{\mathrm{55}}/a_{\mathrm{55}}&0\\ a_{\mathrm{61}}/a_{\mathrm{11}}&a_{\mathrm{62}}/a_{\mathrm{22}}&a_{\mathrm{63}}/a_{\mathrm{33}}&a_{\mathrm{64}}/a_{\mathrm{44}}&a_{\mathrm{65}}/a_{\mathrm{55}}&a_{\mathrm{66}}/a_{\mathrm{66}}\end{array}\right)

Which gives us:

\left(\begin{array}{cccccc}1&0&0&0&0&0\\ a_{\mathrm{21}}/a_{\mathrm{11}}&1&0&0&0&0\\ a_{\mathrm{31}}/a_{\mathrm{11}}&a_{\mathrm{32}}/a_{\mathrm{22}}&1&0&0&0\\ a_{\mathrm{41}}/a_{\mathrm{11}}&a_{\mathrm{42}}/a_{\mathrm{22}}&a_{\mathrm{43}}/a_{\mathrm{33}}&1&0&0\\ a_{\mathrm{51}}/a_{\mathrm{11}}&a_{\mathrm{52}}/a_{\mathrm{22}}&a_{\mathrm{53}}/a_{\mathrm{33}}&a_{\mathrm{54}}/a_{\mathrm{44}}&1&0\\ a_{\mathrm{61}}/a_{\mathrm{11}}&a_{\mathrm{62}}/a_{\mathrm{22}}&a_{\mathrm{63}}/a_{\mathrm{33}}&a_{\mathrm{64}}/a_{\mathrm{44}}&a_{\mathrm{65}}/a_{\mathrm{55}}&1\end{array}\right)

Then replace the first element with -1 and the rest of the elements on the main diagonal with zeros. Call this matrix B:

\left(\begin{array}{cccccc}-1&0&0&0&0&0\\ a_{\mathrm{21}}/a_{\mathrm{11}}&0&0&0&0&0\\ a_{\mathrm{31}}/a_{\mathrm{11}}&a_{\mathrm{32}}/a_{\mathrm{22}}&0&0&0&0\\ a_{\mathrm{41}}/a_{\mathrm{11}}&a_{\mathrm{42}}/a_{\mathrm{22}}&a_{\mathrm{43}}/a_{\mathrm{33}}&0&0&0\\ a_{\mathrm{51}}/a_{\mathrm{11}}&a_{\mathrm{52}}/a_{\mathrm{22}}&a_{\mathrm{53}}/a_{\mathrm{33}}&a_{\mathrm{54}}/a_{\mathrm{44}}&0&0\\ a_{\mathrm{61}}/a_{\mathrm{11}}&a_{\mathrm{62}}/a_{\mathrm{22}}&a_{\mathrm{63}}/a_{\mathrm{33}}&a_{\mathrm{64}}/a_{\mathrm{44}}&a_{\mathrm{65}}/a_{\mathrm{55}}&0\end{array}\right)

Now by matrix multiplication successively square matrix B. That is: Let C=B*B, D=C*C, E=D*D and so on, then the first column in the last matrix (E) will be a convergent. Let those elements be symbolized by the letter r.

\left(\begin{array}{cccccc}r_{\mathrm{11}}&0&0&0&0&0\\ r_{\mathrm{21}}&0&0&0&0&0\\ r_{\mathrm{31}}&0&0&0&0&0\\ r_{\mathrm{41}}&0&0&0&0&0\\ r_{\mathrm{51}}&0&0&0&0&0\\ r_{\mathrm{61}}&0&0&0&0&0\end{array}\right)

Then we only need to divide the ROWS with the elements in the main diagonal of matrix A, to get the first column of the matrix inverse of A.

\left(\begin{array}{cccccc}r_{\mathrm{11}}/a_{\mathrm{11}}&0&0&0&0&0\\ r_{\mathrm{21}}/a_{\mathrm{22}}&0&0&0&0&0\\ r_{\mathrm{31}}/a_{\mathrm{33}}&0&0&0&0&0\\ r_{\mathrm{41}}/a_{\mathrm{44}}&0&0&0&0&0\\ r_{\mathrm{51}}/a_{\mathrm{55}}&0&0&0&0&0\\ r_{\mathrm{61}}/a_{\mathrm{66}}&0&0&0&0&0\end{array}\right)

To calculate the rest of the columns in the inverse of A, repeat for submatrices:

\left(\begin{array}{ccccc}a_{\mathrm{22}}&0&0&0&0\\ a_{\mathrm{32}}&a_{\mathrm{33}}&0&0&0\\ a_{\mathrm{42}}&a_{\mathrm{43}}&a_{\mathrm{44}}&0&0\\ a_{\mathrm{52}}&a_{\mathrm{53}}&a_{\mathrm{54}}&a_{\mathrm{55}}&0\\ a_{\mathrm{62}}&a_{\mathrm{63}}&a_{\mathrm{64}}&a_{\mathrm{65}}&a_{\mathrm{66}}\end{array}\right)

\left(\begin{array}{cccc}a_{\mathrm{33}}&0&0&0\\ a_{\mathrm{43}}&a_{\mathrm{44}}&0&0\\ a_{\mathrm{53}}&a_{\mathrm{54}}&a_{\mathrm{55}}&0\\ a_{\mathrm{63}}&a_{\mathrm{64}}&a_{\mathrm{65}}&a_{\mathrm{66}}\end{array}\right)

and so on.

Edit 21.2.2011: This way of viewing the inverse of a triangular matrix seems to find its explanation when looking at the binomial series expression for inverting a triangular matrix. The inverse of triangular matrix as a binomial series

Mats Granvik mats.granvik(AT)abo.fi

Posted in Uncategorized | 1 Comment

Möbius function identity

x=7; // size of matrix
T=zeros(x,x);

// Value for first element
T(1,1)=1;
T(2,1)=1; // can be any number except 0

// Values for first column
for n=3:x;
T(n,1)=rand(); //random number
end

// Values for second column
for n=2:x;
T(n,2)=T(n,1)-T(n-1,2);
end

// Values for the table from the third column onwards
for n=3:x;
for k=3:n;
s = 0;
for i=1:k-1; 
//recursive definition of divisibility without using the mod function 
s=s+T(n-i,k-1)-T(n-i,k);
end
T(n,k)=s;
end
end
T

A=inv(T)
// Mats Granvik, email: mats.granvik(AT)abo.fi
Posted in Uncategorized | Comments Off on Möbius function identity

Counting primes less than x arithmetically, Scilab, Matlab

// This programs counts the number of primes less than x without using the mod function or any if statements. It could perhaps be said that this is an arithmetic way of counting the primes because here we start with a single one and then sum, subtract, raise to powers, multiply and divide to produce the charachteristic function of the primes.
x=10;
T=zeros(x,x);

// Let the first element be equal to one
T(1,1)=1;
T;
for n=2:x;
for k=2:n;
s = 0;
for i=1:k-1;
//recursive definition of divisibility without using the mod function, table
s=s+T(n-i,k-1)-T(n-i,k);
end
T(n,k)=s;
end
end
T;

A=cumsum(T,'r');

A=inv(A);

B=zeros(x,x);
for n=1:x;
  for k=1:x;
    B(n,k)=(n/k)^(-A(n,k));
  end
end

V=prod(B,'c');

for n=1:x;
  W(n)=V(n)^(-A(n,1));
end

W(1)=0;
for n=2:x;
  W(n)=(W(n)-1)/(n-1);
end

//sequence A000720 in the OEIS, number of primes less than n (or less than x)
pi_n=cumsum(W,'r')
// Mats Granvik, mats.granvik (at) abo.fi
Posted in Uncategorized | 1 Comment

Recursive definition of n mod k, Scilab, Matlab

x=10;  
T=zeros(x,x);  

for n=1:x;
  for k=n+1:x;
    T(n,k)=n;
    end
end

for n=1:x;  
T(n,1)=0;  
end 
T  
for n=2:x;  
for k=2:n;  
s = 0;  
for i=1:k-1;  
//recursive definition of divisibility without using the mod function, table A051731 in the oeis.  
s=s+T(n-i,k-1)-T(n-i,k);  
end 
T(n,k)=s+k-1;  
end 
end 
T  
//Mats Granvik email: mats.granvik(at)abo.fi 
Posted in Uncategorized | Comments Off on Recursive definition of n mod k, Scilab, Matlab