with BitsAKD;
     with ConsAKD; use ConsAKD;
     with ErroAKD;
     with FloaAKD; use FloaAKD;
     with MathAK ;
     with TextAKD; use TextAKD;
     with StriAK;

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

     function  Characteristic ( InFloat: in float ) return long_integer is
     ---------------------------------------------------------------------------
     -- Rounds towards zero.  Example: Characteristic ( -1.999 ) => -1
     -- Long_integer ( float ) will round to the nearest long_integer, or the
     -- nearest _even_ long_integer when it is centered ( 0.5 ==> 0.0 ).
     -- Characteristic ( float ) truncates, instead.
     ---------------------------------------------------------------------------
       Temp : long_integer;
     begin
       if InFloat >= 0.0 then
      Temp := long_integer ( InFloat - 0.5 );
       else
      Temp := long_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;
     exception
       when others =>
      ErroAKD.Notify ( "MathAKD.Characteristic" );
      raise;
     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
     ---------------------------------------------------------------------------
       Temp : float;
     begin
       Temp := MathAK.Factorial ( N )
	 / ( MathAK.Factorial ( R ) * MathAK.Factorial ( N - R ) );
       return Temp;
     exception
       when constraint_error =>
      ErroAKD.Notify ( "MathAKD.Combination", "Constraint_Error" );
      raise;
       when numeric_error =>
      ErroAKD.Notify ( "MathAKD.Combination", "Numeric_Error" );
      raise;
     end Combination;

     procedure Demo is
     ----------------------------------------------------------------------
     begin
       Put_Line ( Copyright );
--     Put_Line ( Description );
       New_Line;
       Put_Line ( "Math sub-programs package" );
       Put_Line ( "for the ""Meridian OpenAda for DOS"" Ada compiler." );
       New_Line;
       Put_Line ( "No demonstration is currently available." );
       New_Line;
       Pause;
     end Demo;

     procedure Float_String (
	      In_Str: in string;
	      Is_Float: out boolean;
	      Out_Float: out float ) is
     ---------------------------------------------------------------------------
       Exponent    : integer := 0;
       Is_Integer  : boolean := false;
       Is_Positive : boolean := true;
       Place       : integer := 0;
       Position    : integer range 1..In_Str'length := 2;
       Temp_Float  : float   := 0.0;
       Temp_Str    : string ( 1..In_Str'length );
     begin
       Is_Float := false;
       Temp_Str := StriAK.Filtered ( In_Str, "_," );
       if Temp_Str ( 1 ) = '-' then
      Is_Positive := false;
       elsif Temp_Str ( 1 ) = '+' then
      Is_Positive := true;
       else
      Is_Positive := true;
      Position := 1;
       end if;
       for index in Position..Temp_Str'last loop
      Position := index;
      exit when Temp_Str ( index ) = '.';
      if Temp_Str ( index ) not in ( '0'..'9' ) then
	return;
      end if;
      Temp_Float := Temp_Float * 10.0
	+ float ( integer'value ( Temp_Str ( index..index ) ) );
       end loop;
       for index in Position + 1 .. Temp_Str'last loop
      Position := index;
      Place := Place + 1;
      exit when Temp_Str ( index ) = 'E' or Temp_Str ( index ) = 'e';
      exit when Temp_Str ( index ) = ASCII.NUL or Temp_Str ( index ) = ' ';
      if Temp_Str ( index ) not in ( '0'..'9' ) then
	return;
      end if;
      Temp_Float := Temp_Float
	+ float ( integer'value ( Temp_Str ( index..index ) ) )
	     / float ( 10 ** Place );
       end loop;
       if Temp_Str ( Position ) = 'E' or Temp_Str ( Position ) = 'e' then
      Integer_String ( Temp_Str ( Position + 1 .. Temp_Str'last ),
	Is_Integer, Exponent );
      if not Is_Integer then
	return;
      end if;
      Temp_Float := Temp_Float * ( 10.0 ** Exponent );
       end if;
       if Is_Positive then
      Out_Float := Temp_Float;
       else
      Out_Float := - Temp_Float;
       end if;
       Is_Float := true;
     exception
       when others =>
      ErroAKD.Notify ( "MathAKD.Float_String ( string; boolean; float )" );
      raise;
     end Float_String;

     procedure Integer_String (
	      In_Str: in string;
	      Is_Integer: out boolean;
	      Out_Integer: out integer ) is
     ---------------------------------------------------------------------------
     begin
       Is_Integer := true;
       Out_Integer := integer'value ( In_Str );
     exception
       when Constraint_Error => Is_Integer := false;
       when others =>
      ErroAKD.Notify ( "MathAKD.Integer_String" );
      raise;
     end Integer_String;

     procedure Integrate
	      ( F: in     TypeAKD.Array_Float_2;
		E: in     float;
		Y: in out TypeAKD.Array_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 TypeAKD;
     ---------------------------------------------------------------------------
     begin
       Y :=  Y + E * ( F * Y );
     exception
       when Constraint_Error =>
      ErroAKD.Notify ( "MathAKD.Integrate", "Constraint_Error" );
      raise;
     end Integrate;

     function  Mantissa ( InFloat: in float ) return float is
     ---------------------------------------------------------------------------
     -- Example:  Mantissa ( 1.3567 ) => 0.3567
     ---------------------------------------------------------------------------
     begin
       return InFloat - float ( MathAKD.Characteristic ( InFloat ) );
     end Mantissa;

     function Mantissa ( InFloat: in float ) return string is
     ---------------------------------------------------------------------------
     -- Example:  Mantissa ( 1.3567 ) ==> ".3567"
     ---------------------------------------------------------------------------
       Mantissa_Tmp: float;
       Mantissa_Str: string ( 1..float'digits + 1 );
     begin
       Mantissa_Tmp := Mantissa ( InFloat );
       Mantissa_Str ( 1 ) := '.';
       for index in 2..Mantissa_Str'last loop
      Mantissa_Tmp := Mantissa ( Mantissa_Tmp ) * 10.0;
      Mantissa_Str ( index ) := StriAK.Image
	( MathAKD.Characteristic ( Mantissa_Tmp ),
	Signed => false ) ( 1 );
       end loop;
       return Mantissa_Str;
     exception
       when others =>
      ErroAKD.Notify ( "MathAKD.Mantissa ( float ) return string" );
      raise;
     end Mantissa;

     function Pos (
       X : in float )
       return float is
     ---------------------------------------------------------------------------
     -- Returns X if it is positive, returns 0.0 otherwise.
     ---------------------------------------------------------------------------
     begin
       if X > 0.0 then
      return X;
       else
      return 0.0;
       end if;
     end Pos;

     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 =>
      ErroAKD.Notify ( "MathAKD.Factorial", "Constraint_Error" );
      raise;
       when numeric_error =>
      ErroAKD.Notify ( "MathAKD.Factorial", "Numeric_Error" );
      raise;
     end Factorial;

     function Sigmoid (
	     InFloat  : in float;
	     Reversed : in boolean := false )
	     return float is
     ---------------------------------------------------------------------------
     --  Input :          - Infinity         0.0        + Infinity
     --  Output:              0.0            0.5            1.0
     ---------------------------------------------------------------------------
     begin
       if Reversed then
      return - Ln ( ( 1.0 / InFloat ) - 1.0 );
       else
      return ( 1.0 / ( 1.0 + ( Exp ( - InFloat ) ) ) );
       end if;
     end Sigmoid;

     function  Sign ( InFloat: in float ) return float is
     ---------------------------------------------------------------------------
     -- Zero will return as a positive.
     ---------------------------------------------------------------------------
     begin
       if InFloat >= 0.0 then
      return +1.0;
       else
      return -1.0;
       end if;
     end Sign;

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

     ---------------------------------------------------------------------------
     ---------------------------------------------------------------------------
     end MathAKD;