with ConsAK ; use ConsAK;
     with FloaAK ; use FloaAK;
     with InteAK ; use InteAK;
     with StriAK ; use StriAK;
     with Text_IO; use Text_IO;
     with TypeAK ; use TypeAK;

     package body MathAK is
     ----------------------------------------------------------------------
     -- Copyright (C) 1994 David Wallace Croft.  All rights reserved.
     ----------------------------------------------------------------------

     function "**" ( X, Y : float )
       return float is
     ----------------------------------------------------------------------
     begin
       return Exp ( Y * Ln ( X ) );
     end "**";

     procedure Alpha (
       V        : in out float;
       dV_dt    : in out float;
       Stimulus : in     float := 0.0;
       Decay    : in     float := 1.0;
       dt       : in     float := 0.001 ) is
     ----------------------------------------------------------------------
     -- Alpha Function                        ___
     -- V(t) = time * exp ( -Decay * time )  /   \______
     -- Time_Peak = 1 / Decay
     ----------------------------------------------------------------------
       dV_dt_2 : float;
     begin
       dV_dt_2 := Stimulus - 2.0 * Decay * dV_dt - ( Decay ** 2 ) * V;
       dV_dt := dV_dt + dt * dV_dt_2;
       V := V + dt * dV_dt;
     end Alpha;

     function Crosscorrelation (
       X, Y : in     Vector_Type )
       return Vector_Type is
     ----------------------------------------------------------------------
     -- page 170 and 479, Introduction to Communications Systems
     -- If X'range = 1..2000, returns -999..1000.
     -- If X'range = 0..1999, returns -1000..999.
     ----------------------------------------------------------------------
       R    : Vector_Type ( X'range ) := ( others => 0.0 );
       Slip : Vector_Type (
	 ( R'last / 2 - R'last + R'first )..( R'last / 2 ) );
       T    : positive;
     begin
       T := X'last - X'first + 1;
       for Tau in R'range loop
	 for Time in X'range loop
	   if Time + Tau <= Y'last then
	     R ( Tau ) := R ( Tau ) + X ( Time ) * Y ( Time + Tau );
	   else
	     R ( Tau )
	       := R ( Tau ) + X ( Time ) * Y ( Time + Tau - Y'last );
	   end if;
	 end loop;
	 R ( Tau ) := R ( Tau ) / float ( T );
       end loop;
       for Index in R'first..( R'last / 2 ) loop
	 Slip ( Index ) := R ( Index );
       end loop;
       for Index in ( R'last / 2 + 1 )..R'last loop
	 Slip ( Index - R'last + R'first - 1 ) := R ( Index );
       end loop;
       return Slip;
     end Crosscorrelation;

     function Characteristic (
       InFloat : in float )
       return integer is
     ----------------------------------------------------------------------
     -- Rounds towards zero.  Example: Characteristic ( -1.999 ) => -1
     -- Integer ( float ) will round to the nearest long_integer, or the
     -- nearest _even_ integer when it is centered ( 0.5 ==> 0.0 ).
     -- Characteristic ( float ) truncates, instead.
     ----------------------------------------------------------------------
       Temp : integer;
     begin
       if ( InFloat >= float ( integer'last ) + 1.0 )
	 or ( InFloat <= float ( integer'first ) - 1.0 ) then
	   raise Constraint_Error;
       end if;
       if Infloat >= 0.0 then
	 Temp := integer ( InFloat - 0.5 );
       else
	 Temp := integer ( InFloat );
       end if;
       if abs ( float ( Temp ) ) + 1.0 = abs ( InFloat ) then
	 if InFloat > 0.0 then
	   Temp := Temp + 1;
	 else
	   Temp := Temp - 1;
	 end if;
       end if;
       return Temp;
     end Characteristic;

     function Combination (
       N : in natural;
       R : in natural )
       return float is
     ---------------------------------------------------------------------------
     -- The number of ways to form combinations of length r from a set of n
     -- distinct objects, repetitions not allowed, is denoted by the symbol
     -- / n \        / n \       n!
     -- |   |, where |   | = ---------------.
     -- \ r /        \ r /   r! * ( n - r )!
     -- "An Introduction to Probability and Its Applications"
     --    by Larsen & Marx, 1985, p117, Theorem 3.3.1.
     -- constraint_error if not N >= R
     ----------------------------------------------------------------------
     begin
       return Factorial ( N ) / ( Factorial ( R ) * Factorial ( N - R ) );
     end Combination;

     function Cos (
       R : in float ) -- radians
       return float is
     ----------------------------------------------------------------------
       Temp  : float := 1.0;
       Temp_New : float;
       R2 : float;
       Index : positive := 2;
       Add   : boolean := false;
       Two_Pi : constant float := 2.0 * Pi;
     begin
       if abs ( R ) > 2.0 * Pi then
	 R2 := R / Two_Pi;
	 R2 := R2 - float ( Characteristic ( R2 ) );
	 R2 := R2 * Two_Pi;
       else
	 R2 := R;
       end if;
       loop
	 begin
	   if Add then
	     Temp_New := Temp + ( R2 ** Index ) / Factorial ( Index );
	   else
	     Temp_New := Temp - ( R2 ** Index ) / Factorial ( Index );
	   end if;
	 exception
	   when others => return Temp;
	 end;
	 exit when Temp = Temp_New;
	 Temp := Temp_New;
	 Index := Index + 2;
	 Add := not Add;
       end loop;
       return Temp;
     end Cos;

     procedure Demo is
     ----------------------------------------------------------------------
       use ConsAK;
       Default :          natural :=  1;
       Min     : constant natural :=  0;
       Max     : constant natural := 13;
       Option  :          natural;
     begin
       loop
	 Put_Line ( Copyright );
	 Put_Line ( "" );
	 Put_Line ( "General math sub-programs package." );
	 Put_Line ( "" );
	 Put_Line ( " 0 = Quit" );
	 Put_Line ( " 1 = ** -- X ** Y, X and Y are of type float" );
	 Put_Line ( " 2 = Characteristic -- truncates float to integer" );
	 Put_Line ( " 3 = Combination" );
	 Put_Line ( " 4 = Cos ( radians )" );
	 Put_Line ( " 5 = Exp -- exponent of an integer" );
	 Put_Line ( " 6 = Exp -- exponent of a float" );
	 Put_Line ( " 7 = Factorial" );
	 Put_Line ( " 8 = Integrate" );
	 Put_Line ( " 9 = Interpolate" );
	 Put_Line ( "10 = Ln -- natural log" );
	 Put_Line ( "11 = Minimum" );
	 Put_Line ( "12 = Odd -- returns true if integer is odd" );
	 Put_Line ( "13 = Sin ( radians )" );
	 Put_Line ( "" );
	 Option := Ask_Nat ( "Option ", Default, Min, Max );
	 Default := Option + 1;
	 if Default > Max then
	   Default := Min;
	 end if;
	 case Option is
	   when  0 =>
	     exit;
	   when  1 =>
	     declare
	       X : float := 4.0;
	       Y : float := 0.5;
	     begin
	       X := Ask ( "X", X );
	       Y := Ask ( "Y", Y );
	       Put_Line ( X ** Y );
	     end;
	     Pause;
	   when  2 =>
	     Put_Line ( Characteristic ( Ask ) );
	     Pause;
	   when  3 =>
	     declare
	       N, R : natural;
	     begin
	       N := Ask_Nat ( "N" );
	       R := Ask_Nat ( "R" );
	       Put_Line ( Combination ( N, R ) );
	     end;
	     Pause;
	   when  4 =>
	     Put_Line ( Cos ( Ask ) );
	     Pause;
	   when  5 =>
	     Put_Line ( Exp ( Ask ) );
	     Pause;
	   when  6 =>
	     Put_Line ( Exp ( Ask ) );
	     Pause;
	   when  7 =>
	     declare
	       F : natural;
	     begin
	       F := Ask_Nat ( "F" );
	       Put_Line ( Factorial ( F ) );
	     end;
	     Pause;
--         when  8 => Integrate;
--         when  9 => Interpolate;
	   when 10 =>
	     Put_Line ( Ln ( Ask ) );
	     Pause;
	   when 11 =>
	     Put_Line ( Minimum ( Ask, Ask ) );
	     Pause;
	   when 12 =>
	     Put_Line ( Odd ( Ask_Int ( "I" ) ) );
	     Pause;
	   when 13 =>
	     Put_Line ( Sin ( Ask ) );
	     Pause;
	   when others =>
	     Put_Line ( "The demonstration of that option is currently"
	       & " not available." );
	     Pause;
	 end case;
	 Put_Line ( "" );
       end loop;
--     exception
--       when others =>
--         ErroAK.Notify ( "MathAK.Demo" );
--         raise;
     end Demo;

     function Exp (
       I : in integer )
       return float is
     ----------------------------------------------------------------------
     begin
       return E ** I;
     end Exp;

     function Exp (
       F : in float )
       return float is
     ----------------------------------------------------------------------
       Index_Max : constant positive := 10;  -- limits accuracy
       Temp   : float := 1.0 + F;
--     Temp_1 : float;
--     Index  : positive := 2;
     begin
       if F = 0.0 then
	 return 1.0;
       elsif F = 1.0 then
	 return E;
       end if;
--     loop
--       begin
--         Temp_1 := Temp + ( F ** Index ) / Factorial ( Index );
--       exception
--         when others => return Temp;
--       end;
--       exit when Temp = Temp_1;
--       Temp := Temp_1;
--       Index := Index + 1;
--     end loop;
       for Index in 2..Index_Max loop
	 Temp := Temp + ( F ** Index ) / Factorial ( Index );
       end loop;
       return Temp;
     exception
       when others =>
	 New_Line;
	 Put_Line ( "Exception in MathAK.Exp!" );
	 New_Line;
	 raise;
     end Exp;

     function Factorial (
       X : in natural )
       return float is
     ----------------------------------------------------------------------
       Result : float := 1.0;
     begin
       for index in 1..X loop
	 Result := Result * float ( index );
       end loop;
       return Result;
       exception
--       when constraint_error =>
--         ErroADK.Notify ( "MathADK.Factorial", "Constraint_Error" );
--         raise;
--       when numeric_error =>
--         ErroADK.Notify ( "MathADK.Factorial", "Numeric_Error" );
--         raise;
	 when others =>
	   New_Line;
	   Put_Line ( "Exception in MathAK.Factorial!" );
	   New_Line;
	   raise;
     end Factorial;

     procedure Integrate (
       Y: in out Float_Array;
       F: in     Float_Array_2;
       E: in     float ) is
     ---------------------------------------------------------------------------
     -- Euler's method
     --
     -- y'' + 0.1y' + y = 0 becomes
     --
     -- y1 = y
     -- y2 = y' =  0.0y1 + 1.0y2     F = ( (  0.0, 1.0 ),
     -- y3 = y''= -1.0y1 - 0.1y2           ( -1.0, -0.1 ) )
     --
     -- Y' = ( y1' , = ( y'   , = FY = F ( y1 , = F ( y   ,
     --        y2' ,     y''  ,            y2 ,       y'  ,
     --        y3' )     y''' )            y3 )       y'' )
     --
     -- Integrate by
     --   Y_New = Y_Old + Epsilon * F * Y_Old
     ---------------------------------------------------------------------------
--     use DC_Type;
     ---------------------------------------------------------------------------
     begin
       null;
--       Y :=  Y + E * ( F * Y );
--     exception
--       when Constraint_Error =>
--         DC_Errr.Notify ( "DC_Math.Integrate", "Constraint_Error" );
--         raise;
     end Integrate;

--   procedure Interpolate (
--     Completed       :    out boolean;
--     Mid_Point       : in out float;
--     Rec             : in out Interpolate_Type;
--     Results         : in     string  := "";
--     Old_File_Name   : in     string  := "";
--     Tmp_File_Name   : in     string  := "" ) is
--   ----------------------------------------------------------------------
--     Completed_Temp : boolean := false;
--   begin
--     if Rec.Start_Over then
--       Open   ( Rec.Old_File, In_File , Old_File_Name );
--       if End_Of_File ( Rec.Old_File ) then
--         Completed_Temp := true;
--         Close ( Rec.Old_File );
--       else
--         Create ( Rec.Tmp_File, Out_File, Tmp_File_Name );
--         Get ( Rec.Old_File, Rec.Old_1 );
--         Get_Line ( Rec.Old_File, Rec.Old_Line_1, Rec.Last_1 );
--       end if;
--     end if;
--     if not Completed_Temp then
--       if not Rec.Start_Over then
--         Put ( Rec.Tmp_File, Mid_Point );
--         Put_Line ( Rec.Tmp_File, Results );
--       else
--         Rec.Start_Over := false;
--       end if;
--       Put ( Rec.Tmp_File, Rec.Old_1 );
--       Put_Line ( Rec.Tmp_File, Rec.Old_Line_1 ( 1..Rec.Last_1 ) );
--       if End_Of_File ( Rec.Old_File ) then
--         Rec.Start_Over := true;
--       else
--         Get ( Rec.Old_File, Rec.Old_2 );
--         Get_Line ( Rec.Old_File, Rec.Old_Line_2, Rec.Last_2 );
--         Mid_Point      := ( Rec.Old_1 + Rec.Old_2 ) / 2.0;
--         Completed_Temp := Mid_Point = Rec.Old_1;
--         Rec.Old_Line_1 := Rec.Old_Line_2;
--         Rec.Last_1     := Rec.Last_2;
--         Rec.Old_1      := Rec.Old_2;
--       end if;
--       if Completed_Temp or Rec.Start_Over then
--         Close ( Rec.Tmp_File );
--         Delete ( Rec.Old_File );
--         Copy ( Tmp_File_Name, Old_File_Name );
--         Delete ( Tmp_File_Name );
--       end if;
--     end if;
--     Completed := Completed_Temp;
--   end Interpolate;
--
--   procedure Interpolate_Demo is
--   ----------------------------------------------------------------------
--     Completed : boolean;
--     Mid_Point : float := 0.0;
--     Rec       : Interpolate_Type;
--     Results   : Line_Type;
--   begin
--     loop
--       Interpolate ( Completed, Mid_Point, Rec,
--         integer'image ( integer ( Mid_Point ) ),
--         "Old_File", "TempFile" );
--       Put ( Mid_Point );
--       Put_Line ( "  " & integer'image ( integer ( Mid_Point ) ) );
--       exit when Completed;
--     end loop;
--   end Interpolate_Demo;

     function Ln (
       X : float )
       return float is
     ----------------------------------------------------------------------
     -- p308, Handbook of Mathematical, Scientific, and Engineering
     --   Formulas, Functions, Tables, Graphs, Transforms
     -- Only valid for X > 0
     ----------------------------------------------------------------------
       Loops : constant integer := 20; -- limits accuracy
       Temp, Temp_Old : float := 0.0;
       Index : positive := 1;
       A : float;
       B : integer;
     begin
       if X < 0.0 then
	 raise Numeric_Error;
       end if;
       A := ( X - 1.0 ) / ( X + 1.0 );
       for Index in 1..Loops loop
	 B := 2 * Index - 1;
	 Temp := Temp + ( A ** B ) / float ( B );
       end loop;
       return 2.0 * Temp;
     end Ln;

     function Minimum (
       X, Y : float )
       return float is
     ----------------------------------------------------------------------
     begin
       if X < Y then
	 return X;
       else
	 return Y;
       end if;
     end Minimum;

     function Odd (
       I : in integer )
       return boolean is
     ----------------------------------------------------------------------
     begin
       return ( I mod 2 ) = 1;
     end Odd;

     function Sigmoid (
       F : float )
       return float is
     ----------------------------------------------------------------------
     begin
       return 1.0 / ( 1.0 + Exp ( -F ) );
     end Sigmoid;

     function Sin (
       R : in float ) -- radians
       return float is
     ----------------------------------------------------------------------
       Temp  : float := R;
       R2    : float;
       Temp_New : float;
       Index : positive := 3;
       Add   : boolean := false;
       Two_Pi : constant float := 2.0 * Pi;
     begin
       if abs ( R ) > 2.0 * Pi then
	 R2 := R / Two_Pi;
	 R2 := R2 - float ( Characteristic ( R2 ) );
	 R2 := R2 * Two_Pi;
	 Temp := R2;
       else
	 R2 := R;
       end if;
       loop
	 begin
	   if Add then
	     Temp_New := Temp + ( R2 ** Index ) / Factorial ( Index );
	   else
	     Temp_New := Temp - ( R2 ** Index ) / Factorial ( Index );
	   end if;
	 exception
	   when others => return Temp;
	 end;
	 exit when Temp = Temp_New;
	 Temp := Temp_New;
	 Index := Index + 2;
	 Add := not Add;
       end loop;
       return Temp;
     end Sin;

     function Sqr (
       F : float )
       return float is
     ----------------------------------------------------------------------
     begin
       return F ** 2;
     end Sqr;

     function  Tanh (
       F :  float )
       return float is
     ---------------------------------------------------------------------------
     --  Input :          - Infinity         0.0        + Infinity
     --  Output:             -1.0            0.0            1.0
     ---------------------------------------------------------------------------
     begin
       return ( Exp ( F ) - Exp ( -F ) ) / ( Exp ( F ) + Exp ( -F ) );
     end Tanh;

     ----------------------------------------------------------------------
     ----------------------------------------------------------------------
     end MathAK;