Return to chebychev_polynomials.adb CVS log

Up to [local] / OpenXM_contrib / PHC / Ada / Schubert

Annotation of OpenXM_contrib/PHC/Ada/Schubert/chebychev_polynomials.adb, Revision 1.1.1.1

1.1       maekawa     1: with Standard_Mathematical_Functions;    use Standard_Mathematical_Functions;
                      2:
                      3: package body Chebychev_Polynomials is
                      4:
                      5:   function Create ( k : natural ) return Vector is
                      6:   begin
                      7:     if k = 0
                      8:      then declare
                      9:             res : Vector(0..0);
                     10:           begin
                     11:             res(0) := 1.0;
                     12:             return res;
                     13:           end;
                     14:      elsif k = 1
                     15:          then declare
                     16:                 res : Vector(0..1);
                     17:               begin
                     18:                 res(0) := 0.0;
                     19:                 res(1) := 1.0;
                     20:                 return res;
                     21:               end;
                     22:          else declare
                     23:                 res : Vector(0..k);
                     24:                 tk1 : constant Vector := Create(k-1);
                     25:                 tk2 : constant Vector := Create(k-2);
                     26:               begin
                     27:                 for i in tk1'range loop
                     28:                   res(i+1) := 2.0*tk1(i);
                     29:                 end loop;
                     30:                 res(0) := 0.0;
                     31:                 for i in tk2'range loop
                     32:                   res(i) := res(i) - tk2(i);
                     33:                 end loop;
                     34:                 return res;
                     35:               end;
                     36:     end if;
                     37:   end Create;
                     38:
                     39:   function Eval ( k : natural; x : double_float ) return double_float is
                     40:
                     41:     res : double_float := COS(double_float(k)*ARCCOS(x));
                     42:
                     43:   begin
                     44:     return res;
                     45:   end Eval;
                     46:
                     47:   function Eval ( p : Vector; x : double_float ) return double_float is
                     48:
                     49:     res : double_float := p(0);
                     50:     pow : double_float := x;
                     51:
                     52:   begin
                     53:     for i in 1..p'last loop
                     54:       res := res + p(i)*pow;
                     55:       pow := pow*x;
                     56:     end loop;
                     57:     return res;
                     58:   end Eval;
                     59:
                     60:   function Diff ( p : Vector ) return Vector is
                     61:
                     62:     res : Vector(0..p'last-1);
                     63:
                     64:   begin
                     65:     for i in res'range loop
                     66:       res(i) := double_float(i+1)*p(i+1);
                     67:     end loop;
                     68:     return res;
                     69:   end Diff;
                     70:
                     71:   function Diff ( p : Vector; k : natural ) return Vector is
                     72:   begin
                     73:     if k = 0
                     74:      then return p;
                     75:      elsif k = 1
                     76:          then return Diff(p);
                     77:          else return Diff(Diff(p),k-1);
                     78:        end if;
                     79:   end Diff;
                     80:
                     81:   function Int ( p : Vector ) return Vector is
                     82:
                     83:     res : Vector(0..p'last+1);
                     84:
                     85:   begin
                     86:     res(0) := 0.0;
                     87:     for i in p'range loop
                     88:       res(i+1) := p(i)/double_float(i+1);
                     89:     end loop;
                     90:     return res;
                     91:   end Int;
                     92:
                     93:   function Int ( p : Vector; k : natural ) return Vector is
                     94:   begin
                     95:     if k = 0
                     96:      then return p;
                     97:      elsif k = 1
                     98:          then return Int(p);
                     99:          else return Int(Int(p),k-1);
                    100:     end if;
                    101:   end Int;
                    102:
                    103: end Chebychev_Polynomials;