// object fourier type 
// a=tlist(['fourier','dim','const','cos','sin'],n,mat(n,1),mat(n,1))
// a_0+a_1*cos(t)+b_1*sin(t)+\cdot+a_m*cos(t)+b_m*sin(t)に対して、
// tlist(['fourier','dim'],m,a_0,
//    (a_1,a_2,\cdots,a_m),(b_1,b_2,\cdots,b_m))を対応させる。

function x=fourier(n)
   s='tlist(["'fourier'","'dim'","'const'","'cos'","'sin'"],n,0,0,0)';
   x=evstr(s);
endfunction


function x=fourier_s(a)
   x=a;
   x('const')=(-1)*a('const');
   x('cos')=(-1)*a('cos');
   x('sin')=(-1)*a('sin');
endfunction


function x=%fourier_a_fourier(a,b)
   c=a('dim');
   d=b('dim');
   if c>=d then
      x=fourier(c);
      x('const')=a('const')+b('const');
      b_cos=mat(c,1);
      b_sin=mat(c,1);
      for i=3:d+2, 
         b_cos(i)=b('cos')(i);
         b_sin(i)=b('sin')(i);
      end; 
      x('cos')=a('cos')+b_cos;
      x('sin')=a('sin')+b_sin;
   else
      x=fourier(d);
      x('const')=a('const')+b('const');
      a_cos=mat(d,1);
      a_sin=mat(d,1);
      for i=3:c+2, 
         a_cos(i)=a('cos')(i);
         a_sin(i)=a('sin')(i);
      end;
      x('cos')=a_cos+b('cos');
      x('sin')=a_sin+b('sin');
   end, 
endfunction


function x=%fourier_s_fourier(a,b)
   c=a('dim');
   d=b('dim');
   if c>=d then
      x=fourier(c);
      x('const')=a('const')-b('const');
      b_cos=mat(c,1);
      b_sin=mat(c,1);
      for i=3:d+2, 
         b_cos(i)=b('cos')(i);
         b_sin(i)=b('sin')(i);
      end; 
      x('cos')=a('cos')-b_cos;
      x('sin')=a('sin')-b_sin;
   else
      x=fourier(d);
      x('const')=a('const')-b('const');
      a_cos=mat(d,1);
      a_sin=mat(d,1);
      for i=3:c+2, 
         a_cos(i)=a('cos')(i);
         a_sin(i)=a('sin')(i);
      end;
      x('cos')=a_cos-b('cos');
      x('sin')=a_sin-b('sin');
   end, 
endfunction


function x=%s_m_fourier(a,b)
   x=b;
   x('const')=a*b('const');
   x('cos')=a*b('cos');
   x('sin')=a*b('sin'); 
endfunction
   

function x=%fourier_m_s(a,b)
   x=a;
   x('const')=b*a('const');
   x('cos')=b*a('cos');
   x('sin')=b*a('sin');
endfunction


function r = fourier_cos_mul(x,m,c)
   s = x('dim');
   s1 = s + m;

   r = fourier(s1);
   r('cos')=mat(s1,1);
   r('sin')=mat(s1,1);
   
   r('cos')(m+2) = r('cos')(m+2) + x('const') * c;
   for i=1:s,
      index = m + i;
      r('cos')(index+2) = r('cos')(index+2) + x('cos')(i+2) * c / 2;
      index = m - i;
      if index == 0 then
          r('const') = r('const') + x('cos')(i+2) * c / 2;
      else
          if index < 0 then
                index = -index;
          end;
          r('cos')(index+2) = r('cos')(index+2) + x('cos')(i+2) * c / 2;
      end,
      
      index = m + i; 
      r('sin')(index+2) = r('sin')(index+2) + x('sin')(i+2) * c / 2;
      index = m - i;
      if index <> 0 then
          if index < 0 then
             index = -index;
             r('sin')(index+2) = r('sin')(index+2) + x('sin')(i+2) * c / 2;
          else
             r('sin')(index+2) = r('sin')(index+2) - x('sin')(i+2) * c / 2;  
          end,
      end,
   end;
endfunction
      

function r = fourier_sin_mul(x,m,c)
   s = x('dim');
   s1 = s + m;

   r = fourier(s1);
   r('cos')=mat(s1,1);
   r('sin')=mat(s1,1);
   
   r('sin')(m+2) = r('sin')(m+2) + x('const') * c;
   for i=1:s,
      index = m + i;
      r('cos')(index+2) = r('cos')(index+2) - x('sin')(i+2) * c / 2;
      index = m - i;
      if index == 0 then
          r('const') = r('const') + x('sin')(i+2) * c / 2;
      else
          if index < 0 then
                index = -index;
          end;
          r('cos')(index+2) = r('cos')(index+2) + x('sin')(i+2) * c / 2;
      end,
      
      index = m + i; 
      r('sin')(index+2) = r('sin')(index+2) + x('cos')(i+2) * c / 2;
      index = m - i;
      if index <> 0 then
          if index < 0 then
             index = -index;
             r('sin')(index+2) = r('sin')(index+2) - x('cos')(i+2) * c / 2;
          else
             r('sin')(index+2) = r('sin')(index+2) + x('cos')(i+2) * c / 2;  
          end,
      end,
   end;
endfunction
   

function r=%fourier_m_fourier(x,y)
   r = x * y('const');
   s = y('dim');

   for i=1:s,
      r1 = fourier_cos_mul(x, i, y('cos')(i+2));
      r = r + r1;
      r2 = fourier_sin_mul(x, i, y('sin')(i+2));
      r = r + r2;
   end,
endfunction


function r=fourier_dif(x)
   r=x;
   s=x('dim');
   r('const')=0;
   for i=1:s,
      r('cos')(i+2) = i * x('sin')(i+2);
      r('sin')(i+2) = (-1)* i * x('cos')(i+2);
   end,
endfunction

function s=fourier_norm(x)
   n=x('dim');
   s=x('const')*x('const');
   for i=1:n,
      s = s + x('cos')(i+2) * x('cos')(i+2)/2;
      s = s + x('sin')(i+2) * x('sin')(i+2)/2;
   end,
   s=sqrt(s);
endfunction    


function t=int_fourier_norm(x)
   n=x('dim');
   s=x('const')*x('const');
   for i=1:n,
      s = s + x('cos')(i+2) * x('cos')(i+2)/2;
      s = s + x('sin')(i+2) * x('sin')(i+2)/2;
   end,
   s = max(abs(s(2)),abs(s(3)))
   t=sqrt(s);
endfunction    


function s=fourier_maxnorm(x)
   n=x('dim');
   s=abs(x('const'));
   for i=1:n,
      s = s + abs(x('cos')(i+2));
      s = s + abs(x('sin')(i+2));
   end,
endfunction    


function s=int_fourier_maxnorm(x)
   n=x('dim');
   s=max(abs(x('const')(2)),abs(x('const')(3)));
   for i=1:n,
      s = s + max(abs(x('cos')(i+2)(2)),abs(x('cos')(i+2)(3)));
      s = s + max(abs(x('sin')(i+2)(2)),abs(x('sin')(i+2)(3)));
   end,
endfunction