program Zagreb;

Uses Croatia;

begin
  Test5;
end.

Unit Croatia;
{
  Cosines of Multiple Angles and Cosine Polynomials
  =================================================
  This program has been designed and is CopyLefted by:

* Han de Bruijn; Systems for Research "A little bit of Physics would be  (===)
* and Education (DTO/SOO), Mekelweg 6  NO Idleness in Mathematics"(HdB) @-O^O-@
* 2628 CD Delft, The Netherlands  http://huizen.dto.tudelft.nl/deBruijn  #/_\#
* E-mail: Han.deBruijn@DTO.TUDelft.NL  Tel: +31 15 27 82751. Fax: 81722   ###

All I ask is to be credited when it is appropriate.
}
INTERFACE

type
  vektor = array of double;
  matrix = array of vektor;

  Cosines = class(TObject)
  private
    function coefficient(L,n : integer) : integer;
    function teken(k : integer) : integer;
    function Inverse(A : matrix) : matrix;
    function Vullen(N : integer; nul : integer) : matrix;
  public
    function Multiple2Powers(N : integer) : matrix;
    function Power2Multiples(N : integer) : matrix;
  end;
{ 
Couple of Tests
}
procedure Test1;
procedure Test2;
procedure serie(CS : Cosines; n : integer);
procedure Test3;
procedure Test4;
procedure Test5;

IMPLEMENTATION

function macht(x : double; n : integer) : double;
{
  Power
}
var
  y : double;
  k : integer;
begin
  y := 1;
  for k := 1 to n do
    y := y * x;
  macht := y;
end;

procedure Test1;
{
  Just a Test
}
const
  N : integer = 33;
var
  k : integer;
  dt,t,u : double;
begin
  dt := 2*Pi/N;
  for k := 0 to N-1 do
  begin
    t := k*dt;
    Writeln(cos(3*t),' = ',4*macht(cos(t),3) - 3*cos(t));
  end;
  Writeln;
  for k := 0 to N-1 do
  begin
    t := k*dt;
    Writeln(cos(4*t),' = ',8*macht(cos(t),4) - 8*macht(cos(t),2) + 1);
  end;
  Writeln;
  for k := 0 to N-1 do
  begin
    t := k*dt;
    u := cos(t);
    Writeln(cos(5*t),' = ',5*u - 20*macht(u,3) + 16*macht(u,5));
  end;
  Writeln;
end;

function boven(n,x : integer) : integer;
{
  C(n,x)
}
var
  k,teller,noemer : integer;
begin
  teller := 1;
  noemer := 1;
  for k := 1 to x do
  begin
    teller := teller * (n-k+1);
    noemer := noemer * k;
    if (teller mod noemer) = 0 then
    begin
      teller := (teller div noemer);
      noemer := 1;
    end;
  end;
  boven := teller;
end;

function Cosines.coefficient(L,n : integer) : integer;
{
  of Cosine Polynomial
}
var
  m,som : integer;
begin
  som := 0;
  for m := L to (n div 2) do
  begin
    som := som + boven(n,2*m) * boven(m,L);
  end;
  coefficient := som;
end;

procedure Test2;
{
  Just a Test
}
var
  k : integer;
  CS : Cosines;
begin
  CS := Cosines.Create;
  Write(CS.coefficient(0,1):4);
  Writeln;
  Write(CS.coefficient(0,2):4);
  Write(CS.coefficient(1,2):4);
  Writeln;
  Write(CS.coefficient(0,3):4);
  Write(CS.coefficient(1,3):4);
  Writeln;
  for k := 0 to 2 do
    Write(CS.coefficient(k,4):4);
  Writeln;
  for k := 0 to 2 do
    Write(CS.coefficient(k,5):4);
  Writeln;
end;

function Cosines.teken(k : integer) : integer;
{
  (-1)^k
}
var
  t : integer;
begin
  t := -1;
  if (k mod 2) = 0 then t := +1;
  teken := t;
end;

procedure serie(CS : Cosines; n : integer);
{
  Cosine Polynomial
}
var
  L,t,m : integer;
begin
  for L := (n div 2) downto 0 do
  begin
    t := CS.teken(L);
    if t < 0 then Write(' - ');
    if t > 0 then Write(' + ');
    Write(CS.coefficient(L,n));
    m := n-2*L;
    if m = 1 then Write('.cos(x)');    
    if m > 1 then Write('.cos^',m,'(x)');    
  end;
  Writeln;
end;

procedure Test3;
{
  Just a Test
}
var
  k : integer;
  CS : Cosines;
begin
  CS := Cosines.Create;
  for k := 1 to 9 do
    serie(CS,k);
end;

function Cosines.Inverse(A : matrix) : matrix;
{
  Matrix Inversion
}
var
  b : matrix;
  pe : double;
  ndm,i,j,k : integer;
begin
  ndm := Length(A);
  SetLength(b,ndm,ndm);
  for i := 0 to ndm-1 do
    for j := 0 to ndm-1 do
      b[i,j] := A[i,j];

  for k := 0 to ndm-2 do
  begin
    pe := b[k,k];
    for i := k+1 to ndm-1 do
    begin
      b[i,k] := - b[i,k]/pe;
      for j := 0 to k-1 do
        b[i,j] := b[i,j] + b[i,k]*b[k,j];
    end;
  end;
  for i := 0 to ndm-1 do
  begin
    pe := b[i,i];
    b[i,i] := 1/pe;
    for j := 0 to i-1 do
      b[i,j] := b[i,j]/pe;
  end;
  Inverse := b;
end;

procedure Afdrukken(h : matrix);
{
  Matrix Print
}
var
  ndm,i,j : integer;
begin
  ndm := Length(h);
  for i := 0 to ndm-1 do
  begin
    for j := 0 to ndm-1 do
    begin
      Write(h[i][j]:9:3);
    end;
    Writeln;
  end;
end;

function Cosines.Vullen(N : integer; nul : integer) : matrix;
{
  Fill Matrices
}
var
  even : boolean;
  ndm,i,j : integer;
  h : matrix;
begin
  ndm := ((N+1) div 2);
  even := ((N mod 2) = 0) and (nul = 0);
  if even then ndm := ndm + 1;
  SetLength(h,ndm,ndm);
  for i := 0 to ndm-1 do
  for j := 0 to ndm-1 do
    h[i][j] := 0;
    
  for i := 0 to ndm-1 do
  begin
    for j := i downto 0 do
      h[i][i-j] := teken(j)*coefficient(j,2*i+nul);
  end;
  Vullen := h;
end;

function Maal(A,B : matrix) : matrix;
{
  Matrix Multiply
}
var
  ndm,i,j,k : integer;
  C : matrix;
  h : double;
begin
  ndm := Length(A);
  SetLength(C,ndm,ndm);
  for i := 0 to ndm-1 do
  begin
    for j := 0 to ndm-1 do
    begin
      h := 0;
      for k := 0 to ndm-1 do
        h := h + A[i,k] * B[k,j];
      C[i,j] := h;
    end;
  end;
  Maal := C;
end;
  
procedure Test4;
{
  Just a Test
}
var
  A,B,C : matrix;
  CS : Cosines;
begin
  CS := Cosines.Create;
  SetLength(A,0,0);
  SetLength(B,0,0);
  SetLength(C,0,0);
  A := CS.Vullen(7,0);
  Afdrukken(A);
  Writeln;
  B := CS.Inverse(A);
  Afdrukken(B);
  Writeln;
  C := Maal(A,B);
  Afdrukken(C);
  Writeln;
  A := CS.Vullen(7,1);
  Afdrukken(A);
  Writeln;
  B := CS.Inverse(A);
  Afdrukken(B);
  Writeln;
  C := Maal(A,B);
  Afdrukken(C);
end;

function Cosines.Multiple2Powers(N : integer) : matrix;
{
  Expand Cosine of Multiple Angle
  as a Polynomial of Cosine Powers 
}
var
  i,j,M : integer;
  A,C : matrix;
begin
  SetLength(C,N+1,N+1);
  for i := 0 to N do
    for j := 0 to N do
      C[i,j] := 0; 

{ Even 0 .. N }
  A := Vullen(N,0);
  M := Length(A);
  for i := 0 to M-1 do
  begin
    for j := 0 to i do
    begin
      C[2*i,2*j] := A[i,j];     
    end;
  end;
{ Odd 1 .. N }
  A := Vullen(N,1);
  M := Length(A);
  for i := 0 to M-1 do
  begin
    for j := 0 to i do
    begin
      C[2*i+1,2*j+1] := A[i,j];     
    end;
  end;
  Multiple2Powers := C;
end;

function Cosines.Power2Multiples(N : integer) : matrix;
{
  Expand Cosine Power as a Series
  of Cosines of a Multiple Angle
}
var
  i,j,M : integer;
  A,C : matrix;
begin
  SetLength(C,N+1,N+1);
  for i := 0 to N do
    for j := 0 to N do
      C[i,j] := 0; 

{ Even 0 .. N }
  A := Vullen(N,0);
  M := Length(A);
  A := Inverse(A);
  for i := 0 to M-1 do
  begin
    for j := 0 to i do
    begin
      C[2*i,2*j] := A[i,j];     
    end;
  end;
{ Odd 0 .. N }
  A := Vullen(N,1);
  M := Length(A);
  A := Inverse(A);
  for i := 0 to M-1 do
  begin
    for j := 0 to i do
    begin
      C[2*i+1,2*j+1] := A[i,j];     
    end;
  end;
  Power2Multiples := C;
end;

procedure Test5;
{
  Just a Test
}
var
  A : matrix;
  CS : Cosines;
begin
  CS := Cosines.Create;
  A := CS.Multiple2Powers(7); 
  Afdrukken(A);
  Writeln;
  A := CS.Power2Multiples(7); 
  Afdrukken(A);
end;

END.