package body GENERIC_DIGITS_ARITHMETIC is

  BASE : constant := 10 ;
  HALF_BASE : constant := BASE / 2 ;
  CONSTANT_ZERO : D_NUMBER ;
  CONSTANT_ONE : D_NUMBER ;
  CONSTANT_TWO : D_NUMBER ;
  CONSTANT_THREE : D_NUMBER ;
  CONSTANT_FOUR : D_NUMBER ;
  CONSTANT_SIX : D_NUMBER ;
  CONSTANT_EIGHT : D_NUMBER ;
  CONSTANT_TEN : D_NUMBER ;
  CONSTANT_TWENTY_FOUR : D_NUMBER ;
  CONSTANT_EPSILON : D_NUMBER ;
  CONSTANT_HALF : D_NUMBER ;
  CONSTANT_FOURTH : D_NUMBER ;
  CONSTANT_EIGHTH : D_NUMBER ;
  CONSTANT_PI : D_NUMBER ;
  CONSTANT_TWO_PI : D_NUMBER ;
  CONSTANT_PI_OVER_TWO : D_NUMBER  ;
  CONSTANT_PI_OVER_THREE : D_NUMBER  ;
  CONSTANT_PI_OVER_FOUR : D_NUMBER  ;
  CONSTANT_PI_OVER_SIX : D_NUMBER  ;
  CONSTANT_NATURAL_E : D_NUMBER  ;
  LOG_BASE_E_OF_10 : D_NUMBER ;
  LOG_BASE_E_OF_2 : D_NUMBER ;
  CONSTANT_SQRT_3 : D_NUMBER ;
  CONSTANT_ONE_OVER_SQRT_3 : D_NUMBER ;
  CONSTANT_ATAN_BREAK_1 : D_NUMBER ;
  SQRT_ATAN_BREAK_1 : D_NUMBER ;
  CONSTANT_ATAN_BREAK_2 : D_NUMBER ;
  FACTORIAL : array ( 0 .. INTEGER(SIZE'LAST) ) of D_NUMBER ;
  ATAN_C : array ( 0 .. (INTEGER(SIZE'LAST)*4)/3) of D_NUMBER ;
  ONE_OVER : array ( 0 .. (INTEGER(SIZE'LAST)*5)/2 ) of D_NUMBER ;

  package INT_IO is new TEXT_IO.INTEGER_IO ( INTEGER ) ;

  function DIGITS_ZERO return D_NUMBER is
    ZERO : D_NUMBER ;  -- zero by default, supplied for consistancy
  begin
    return ZERO ;
  end DIGITS_ZERO ;

  function DIGITS_ONE return D_NUMBER is
    ONE : D_NUMBER ;
  begin
    ONE.EXP := 1 ;
    ONE.D_DIGITS ( SIZE'LAST ) := 1 ;
    return ONE ; -- basically 0.1 * 10.0 ** 1 = 1.0
  end DIGITS_ONE ;

  function DIGITS_EPSILON return D_NUMBER is
  begin
    return CONSTANT_EPSILON ; -- basically 0.1 * 10.0 ** (-SIZE'LAST+1)
  end DIGITS_EPSILON ;

  function DIGITS_LAST return D_NUMBER is
    LAST : D_NUMBER := CONSTANT_ONE - DIGITS_EPSILON ;
  begin
    LAST.EXP := INTEGER'LAST ;
    return LAST ; -- basically 0.999... * 10.0 ** INTEGER'LAST-1
  end DIGITS_LAST ;

  function DIGITS_FIRST return D_NUMBER is
  begin
    return - DIGITS_LAST ;
  end DIGITS_FIRST ;

  function DIGITS_SMALL return D_NUMBER is
  begin
    return CONSTANT_ONE / DIGITS_LAST ; -- about 10.0 ** (-INTEGER'LAST+1)
  end DIGITS_SMALL ;

  function "+" ( LEFT , RIGHT : D_NUMBER ) return D_NUMBER is
    SUM : D_NUMBER ;
    CARRY : INTEGER := 0 ;
    OFFSET : INTEGER ;
  begin
    if LEFT.D_DIGITS ( SIZE'LAST ) = 0 then
      return RIGHT ;
    elsif RIGHT.D_DIGITS ( SIZE'LAST ) = 0 then
      return LEFT ;
    elsif RIGHT.EXP - LEFT.EXP > INTEGER ( SIZE'LAST ) + 1 then
      return RIGHT ;
    elsif LEFT.EXP - RIGHT.EXP > INTEGER ( SIZE'LAST ) + 1 then
      return LEFT ;
    end if ;
    if LEFT.SIGN = RIGHT.SIGN then
      if LEFT.EXP >= RIGHT.EXP then
        SUM := LEFT ;
        OFFSET := LEFT.EXP - RIGHT.EXP ;
        if OFFSET > 0 and then
           OFFSET <= INTEGER ( SIZE'LAST ) + 1 and then
           RIGHT.D_DIGITS ( SIZE( OFFSET - 1 )) >= HALF_BASE then
           CARRY := 1 ;  -- round
        end if ;
        for I in SIZE loop
          if INTEGER ( I ) + OFFSET <= INTEGER ( SIZE'LAST ) then
            SUM.D_DIGITS ( I ) := SUM.D_DIGITS ( I ) + RIGHT.D_DIGITS ( I +
               SIZE( OFFSET )) ;
          end if ;
          SUM.D_DIGITS ( I ) := SUM.D_DIGITS ( I ) + CARRY ;
          if SUM.D_DIGITS ( I ) >= BASE then
            SUM.D_DIGITS ( I ) := SUM.D_DIGITS ( I ) - BASE ;
            CARRY := 1 ;
          else
            CARRY := 0 ;
          end if ;
        end loop ;
      else
        SUM := RIGHT ;
        OFFSET := RIGHT.EXP - LEFT.EXP ;
        if OFFSET > 0 and then
           OFFSET <= INTEGER ( SIZE'LAST ) + 1 and then
           LEFT.D_DIGITS ( SIZE( OFFSET - 1 )) >= HALF_BASE then
          CARRY := 1 ;  -- round
        end if ;
        for I in SIZE loop
          if INTEGER ( I ) + OFFSET <= INTEGER ( SIZE'LAST ) then
            SUM.D_DIGITS ( I ) := SUM.D_DIGITS ( I ) + LEFT.D_DIGITS ( I +
               SIZE( OFFSET )) ;
          end if ;
          SUM.D_DIGITS ( I ) := SUM.D_DIGITS ( I ) + CARRY ;
          if SUM.D_DIGITS ( I ) >= BASE then
            SUM.D_DIGITS ( I ) := SUM.D_DIGITS ( I ) - BASE ;
            CARRY := 1 ;
          else
            CARRY := 0 ;
          end if ;
        end loop ;
      end if ;
    else  -- signs differ, use "-" to do operation
      if LEFT.SIGN = PLUS then
        SUM := RIGHT ;
        SUM.SIGN := PLUS ;
        return LEFT - SUM ;
      else
        SUM := LEFT ;
        SUM.SIGN := PLUS ;
        return RIGHT - SUM ;
      end if ;
    end if ;

-- clean up if carry out of top digit
    if CARRY = 1 then
      SUM.EXP := SUM.EXP + 1 ;
      if SUM.D_DIGITS ( 0 ) >= HALF_BASE then  -- round
        for I in SIZE ( 1 ) .. SIZE'LAST loop
          SUM.D_DIGITS ( I - 1 ) := SUM.D_DIGITS ( I ) + CARRY ;
          if SUM.D_DIGITS ( I - 1 ) >= BASE then
            SUM.D_DIGITS ( I - 1 ) := SUM.D_DIGITS ( I - 1 ) - BASE ;
            CARRY := 1 ;
          else
            CARRY := 0 ;
          end if ;
        end loop ;
        SUM.D_DIGITS ( SIZE'LAST ) := CARRY + 1 ;
      else
        for I in SIZE ( 1 ) .. SIZE'LAST loop
          SUM.D_DIGITS ( I - 1 ) := SUM.D_DIGITS ( I ) ;
        end loop ;
        SUM.D_DIGITS ( SIZE'LAST ) := 1 ;
      end if ;
    end if ;
    return SUM ;
  end "+" ;

  function "-" ( LEFT , RIGHT : D_NUMBER ) return D_NUMBER is
    DIFF : D_NUMBER ;
    BORROW : INTEGER := 0 ;
    OFFSET : INTEGER ;
  begin
    if LEFT.D_DIGITS ( SIZE'LAST ) = 0 and RIGHT.D_DIGITS(SIZE'LAST) /= 0 then
      DIFF := RIGHT ;
      if RIGHT.SIGN = PLUS then
        DIFF.SIGN := MINUS ;
      else
        DIFF.SIGN := PLUS ;
      end if ;
      return DIFF ;
    elsif RIGHT.D_DIGITS ( SIZE'LAST ) = 0 then
      return LEFT ;
    elsif RIGHT.EXP - LEFT.EXP > INTEGER ( SIZE'LAST ) + 1 then
      return RIGHT ;
    elsif LEFT.EXP - RIGHT.EXP > INTEGER ( SIZE'LAST ) + 1 then
      return LEFT ;
    end if ;
    if LEFT.SIGN = RIGHT.SIGN then
      if abs LEFT >= abs RIGHT then
        DIFF := LEFT ;
        OFFSET := LEFT.EXP - RIGHT.EXP ;
        if OFFSET > 0 and then
           OFFSET <= INTEGER ( SIZE'LAST ) + 1 and then
           RIGHT.D_DIGITS ( SIZE( OFFSET - 1 )) >= HALF_BASE then
          BORROW := 1 ;  -- round
        end if ;
        for I in SIZE loop
          if INTEGER ( I ) + OFFSET <= INTEGER ( SIZE'LAST ) then
            DIFF.D_DIGITS ( I ) := DIFF.D_DIGITS ( I ) - RIGHT.D_DIGITS ( I +
               SIZE( OFFSET )) ;
          end if ;
          DIFF.D_DIGITS ( I ) := DIFF.D_DIGITS ( I ) - BORROW ;
          if DIFF.D_DIGITS ( I ) < 0 then
            DIFF.D_DIGITS ( I ) := DIFF.D_DIGITS ( I ) + BASE ;
            BORROW := 1 ;
          else
            BORROW := 0 ;
          end if ;
        end loop ;
      else
        DIFF := RIGHT ;
        if RIGHT.SIGN = PLUS then
          DIFF.SIGN := MINUS ;
        else
          DIFF.SIGN := PLUS ;
        end if ;
        OFFSET := RIGHT.EXP - LEFT.EXP ;
        if OFFSET > 0 and then
           OFFSET <= INTEGER ( SIZE'LAST ) + 1 and then
           LEFT.D_DIGITS ( SIZE( OFFSET - 1 )) >= HALF_BASE then
          BORROW := 1 ;  -- round
        end if ;
        for I in SIZE loop
          if INTEGER ( I ) + OFFSET <= INTEGER ( SIZE'LAST ) then
            DIFF.D_DIGITS ( I ) := DIFF.D_DIGITS ( I ) - LEFT.D_DIGITS ( I +
               SIZE( OFFSET )) ;
          end if ;
          DIFF.D_DIGITS ( I ) := DIFF.D_DIGITS ( I ) - BORROW ;
          if DIFF.D_DIGITS ( I ) < 0 then
            DIFF.D_DIGITS ( I ) := DIFF.D_DIGITS ( I ) + BASE ;
            BORROW := 1 ;
          else
            BORROW := 0 ;
          end if ;
        end loop ;
      end if ;
    else  -- signs differ, use "+" to do operation
      if LEFT.SIGN = PLUS then
        DIFF := RIGHT ;
        DIFF.SIGN := PLUS ;
        return LEFT + DIFF ;
      else
        DIFF := RIGHT ;
        DIFF.SIGN := MINUS ;
        return LEFT + DIFF ;
      end if ;
    end if ;

-- clean up if top digit is zero
    if DIFF.D_DIGITS ( SIZE'LAST ) = 0 then
      DIFF := NORMALIZE ( DIFF ) ;
    end if ;
    return DIFF ;
  end "-" ;

  function "*" ( LEFT , RIGHT : D_NUMBER ) return D_NUMBER is
    PROD : D_NUMBER ;
    N : INTEGER := 2 * INTEGER ( SIZE'LAST ) + 1 ;
    SUM : array ( 0 .. N ) of INTEGER ;
    TEST : INTEGER ;
    K : INTEGER ;
  begin
    if LEFT.D_DIGITS ( SIZE'LAST ) = 0 or
       RIGHT.D_DIGITS ( SIZE'LAST ) = 0 then
      return PROD ;
    end if ;
    if LEFT.EXP < 0 and
       RIGHT.EXP < 0 then
      begin
        TEST := LEFT.EXP + RIGHT.EXP ;
        TEST := TEST - 1 ;
      exception
        when CONSTRAINT_ERROR | NUMERIC_ERROR =>
          return PROD ;
      end ;  -- underflow check block
      if TEST = INTEGER'FIRST then
        return PROD ;
      end if ;
    end if ;
    for I in 0 .. N loop
      SUM ( I ) := 0 ;
    end loop ;
    for I in SIZE loop
      for J in SIZE loop
        K := INTEGER ( I ) + INTEGER ( J ) ;
        SUM ( K ) := SUM ( K ) + LEFT.D_DIGITS ( J ) * RIGHT.D_DIGITS ( I ) ;
        if SUM ( K ) >= BASE then
          SUM ( K + 1 ) := SUM ( K + 1 ) + SUM ( K ) / BASE ;
          SUM ( K ) := SUM ( K ) mod BASE ;
        end if ;
      end loop ;
    end loop ;
    if LEFT.SIGN /= RIGHT.SIGN then
      PROD.SIGN := MINUS ;
    end if ;
    if SUM ( N ) /= 0 then
      PROD.EXP := LEFT.EXP + RIGHT.EXP ;
      for I in SIZE loop
        PROD.D_DIGITS ( I ) := SUM ( INTEGER( SIZE'LAST ) + 1 + INTEGER ( I )
          ) ;
      end loop ;
    else
      PROD.EXP := LEFT.EXP + RIGHT.EXP - 1 ;
      for I in SIZE loop
        PROD.D_DIGITS ( I ) := SUM ( INTEGER( SIZE'LAST ) + INTEGER ( I )) ;
      end loop ;
    end if ;
    return PROD ;
  end "*" ;

  function "/" ( LEFT , RIGHT : D_NUMBER ) return D_NUMBER is
    QUO : D_NUMBER ;
    N : INTEGER := INTEGER ( SIZE'LAST ) ;
    QUOTIENT : array ( INTEGER( - 1 ) .. N) of INTEGER ;
    DIVISOR : array ( 0 .. N ) of INTEGER ;
    REMAINDER : array ( 0 .. 2 * N + 1 ) of INTEGER ;
    TEST : INTEGER ;
    REMAINDER_INDEX : INTEGER := 2 * N + 1 ;
    REMAINDER_TRIAL : INTEGER ;
  begin
    if RIGHT.D_DIGITS ( SIZE'LAST ) = 0 then
      raise NUMERIC_ERROR ;
    elsif LEFT.D_DIGITS ( SIZE'LAST ) = 0 then
      return QUO ;
    end if ;
    if LEFT.EXP < 0 and
       RIGHT.EXP > 0 then
      begin
        TEST := LEFT.EXP - RIGHT.EXP ;
        TEST := TEST - 1 ;
      exception
        when CONSTRAINT_ERROR | NUMERIC_ERROR =>
          return QUO ;
      end ;  -- underflow check block
      if TEST = INTEGER'FIRST then
        return QUO ;
      end if ;
    end if ;
    for I in 0 .. N loop
      DIVISOR ( I ) := RIGHT.D_DIGITS ( SIZE( I )) ;
      REMAINDER ( I ) := 0 ;
      REMAINDER ( I + N + 1 ) := LEFT.D_DIGITS ( SIZE( I )) ;
    end loop ;
    for L in reverse INTEGER ( - 1 ) .. N loop
      REMAINDER_TRIAL := REMAINDER ( REMAINDER_INDEX ) ;
      if REMAINDER_INDEX < 2 * N + 1 then
        REMAINDER_TRIAL := REMAINDER_TRIAL + BASE * REMAINDER (
           REMAINDER_INDEX + 1 ) ;
      end if ;
      QUOTIENT ( L ) := REMAINDER_TRIAL / DIVISOR ( N ) ;
      if QUOTIENT ( L ) > 0 then
        for I in 0 .. N loop
          REMAINDER ( L + I + 1 ) := REMAINDER ( L + I + 1 ) - QUOTIENT ( L )
             * DIVISOR ( I ) ;
          while REMAINDER ( L + I + 1 ) < 0 and L + I + 2 <= 2 * N + 1 loop
            REMAINDER ( L + I + 2 ) := REMAINDER ( L + I + 2 ) - 1 ;
            REMAINDER ( L + I + 1 ) := REMAINDER ( L + I + 1 ) + BASE ;
          end loop ;
        end loop ;
        if L = N then
          while REMAINDER ( L + N + 1 ) < 0 loop  -- add in divisor
            QUOTIENT ( L ) := QUOTIENT ( L ) - 1 ;
            for I in 0 .. N loop
              REMAINDER ( L + I + 1 ) := REMAINDER ( L + I + 1 ) + DIVISOR ( I
                 ) ;
              if REMAINDER ( L + I + 1 ) >= BASE and
                 L + I + 2 <= 2 * N + 1 then
                REMAINDER ( L + I + 2 ) := REMAINDER ( L + I + 2 ) + 1 ;
                REMAINDER ( L + I + 1 ) := REMAINDER ( L + I + 1 ) - BASE ;
              end if ;
            end loop ;
          end loop ;
        else
          while REMAINDER ( L + N + 2 ) < 0 loop  -- add in divisor
            QUOTIENT ( L ) := QUOTIENT ( L ) - 1 ;
            for I in 0 .. N loop
              REMAINDER ( L + I + 1 ) := REMAINDER ( L + I + 1 ) + DIVISOR ( I
                 ) ;
              if REMAINDER ( L + I + 1 ) >= BASE and
                 L + I + 2 <= 2 * N + 1 then
                REMAINDER ( L + I + 2 ) := REMAINDER ( L + I + 2 ) + 1 ;
                REMAINDER ( L + I + 1 ) := REMAINDER ( L + I + 1 ) - BASE ;
              end if ;
            end loop ;
          end loop ;
        end if ;
      end if ;
      REMAINDER_INDEX := REMAINDER_INDEX - 1 ;
    end loop ;
    if LEFT.SIGN /= RIGHT.SIGN then
      QUO.SIGN := MINUS ;
    end if ;
    QUO.EXP := LEFT.EXP - RIGHT.EXP ;
    if QUOTIENT ( N ) /= 0 then
      QUO.EXP := QUO.EXP + 1 ;
      for I in SIZE loop
        QUO.D_DIGITS ( I ) := QUOTIENT ( INTEGER( I )) ;
      end loop ;
    else
      for I in SIZE loop
        QUO.D_DIGITS ( I ) := QUOTIENT ( INTEGER( I ) - 1) ;
      end loop ;
    end if ;
    return QUO ;
  end "/" ;

  function "**" ( LEFT : D_NUMBER ;
                  RIGHT : INTEGER ) return D_NUMBER is
    POWER : D_NUMBER := LEFT ;
    COUNT : INTEGER := abs RIGHT - 1 ;
  begin
    if RIGHT = 0 then
      return CONSTANT_ONE ;
    elsif RIGHT = 1 then
      return LEFT ;
    elsif RIGHT < 0 and LEFT.D_DIGITS ( SIZE'LAST ) = 0 then
      raise NUMERIC_ERROR ;
    elsif RIGHT = -1 then
      return CONSTANT_ONE / LEFT ;
    end if ;
    for I in 1..COUNT loop
      POWER := POWER * LEFT ;
    end loop ;
    if RIGHT < 0 then
      POWER := CONSTANT_ONE / POWER ;
    end if ;
    return POWER ;
  end "**" ;

  function "+" ( RIGHT : D_NUMBER ) return D_NUMBER is
  begin
    return RIGHT ;
  end "+" ;

  function "-" ( RIGHT : D_NUMBER ) return D_NUMBER is
    RESULT : D_NUMBER := RIGHT ;
  begin
    if RESULT.D_DIGITS ( SIZE'LAST ) = 0 then
      return RESULT ;
    end if ;
    if RESULT.SIGN = PLUS then
      RESULT.SIGN := MINUS ;
    else
      RESULT.SIGN := PLUS ;
    end if ;
    return RESULT ;
  end "-" ;

  function "abs" ( RIGHT : D_NUMBER ) return D_NUMBER is
    RESULT : D_NUMBER := RIGHT ;
  begin
    RESULT.SIGN := PLUS ;
    return RESULT ;
  end "abs" ;

  function ">" ( LEFT , RIGHT : D_NUMBER ) return BOOLEAN is
  begin
    if LEFT = RIGHT then
      return FALSE ;
    end if ;
    if LEFT.SIGN = MINUS and
       RIGHT.SIGN = PLUS then
      return FALSE ;
    elsif LEFT.SIGN = PLUS and
       RIGHT.SIGN = MINUS then
      return TRUE ;
    elsif LEFT.D_DIGITS ( SIZE'LAST ) = 0 and
       RIGHT.SIGN = PLUS then
      return FALSE ;
    elsif RIGHT.D_DIGITS ( SIZE'LAST ) = 0 and
       LEFT.SIGN = PLUS then
      return TRUE ;
    end if ;
    if LEFT.SIGN = PLUS then  -- RIGHT.SIGN = PLUS
      if LEFT.EXP > RIGHT.EXP then
        return TRUE ;
      elsif RIGHT.EXP > LEFT.EXP then
        return FALSE ;
      else  -- exponents equal, check digits
        for I in reverse SIZE loop
          if LEFT.D_DIGITS ( I ) > RIGHT.D_DIGITS ( I ) then
            return TRUE ;
          elsif LEFT.D_DIGITS ( I ) < RIGHT.D_DIGITS ( I ) then
            return FALSE ;
          end if ;
        end loop ;
        return FALSE ;
      end if ;
    else  -- LEFT.SIGN = MINUS and RIGHT.SIGN = MINUS
      if LEFT.EXP > RIGHT.EXP then
        return FALSE ;
      elsif RIGHT.EXP > LEFT.EXP then
        return TRUE ;
      else  -- exponents equal, check digits
        for I in reverse SIZE loop
          if LEFT.D_DIGITS ( I ) > RIGHT.D_DIGITS ( I ) then
            return FALSE ;
          elsif LEFT.D_DIGITS ( I ) < RIGHT.D_DIGITS ( I ) then
            return TRUE ;
          end if ;
        end loop ;
        return FALSE ;
      end if ;
    end if ;
  end ">" ;

  function ">=" ( LEFT , RIGHT : D_NUMBER ) return BOOLEAN is
  begin
    if LEFT = RIGHT then
      return TRUE ;
    end if ;
    return LEFT > RIGHT ;
  end ">=" ;

  function "<" ( LEFT , RIGHT : D_NUMBER ) return BOOLEAN is
  begin
    return RIGHT > LEFT ;
  end "<" ;

  function "<=" ( LEFT , RIGHT : D_NUMBER ) return BOOLEAN is
  begin
    if LEFT = RIGHT then
      return TRUE ;
    end if ;
    return RIGHT > LEFT ;
  end "<=" ;

  function D_NUMBER_TO_INTEGER ( ITEM : D_NUMBER ) return INTEGER is
    NUMBER : INTEGER := 0 ;
    POWER : INTEGER := BASE ;
  begin
    if ITEM.EXP <= 0 then
      return 0 ;
    end if ;
    NUMBER := ITEM.D_DIGITS ( SIZE'LAST - SIZE( ITEM.EXP - 1 )) ;
    if ITEM.EXP > 1 then
      for I in SIZE'LAST - SIZE ( ITEM.EXP - 2 ) .. SIZE'LAST loop
        NUMBER := NUMBER + POWER * ITEM.D_DIGITS ( I ) ;
        POWER := POWER * BASE ;
      end loop ;
    end if ;
    if ITEM.SIGN = MINUS then
      NUMBER := - NUMBER ;
    end if ;
    return NUMBER ;
  end D_NUMBER_TO_INTEGER ;

  function INTEGER_TO_D_NUMBER ( ITEM : INTEGER ) return D_NUMBER is
    NUMBER : D_NUMBER ;
    REST : INTEGER ;
  begin
    if ITEM < 0 then
      NUMBER.SIGN := MINUS ;
    end if ;
    REST := abs ITEM ;
    NUMBER.EXP := INTEGER ( SIZE'LAST ) + 1 ;
    for I in SIZE loop
      exit when REST = 0 ;
      NUMBER.D_DIGITS ( I ) := REST mod BASE ;
      REST := REST / BASE ;
    end loop ;
    return NORMALIZE ( NUMBER ) ;
  end INTEGER_TO_D_NUMBER ;

  function D_NUMBER_TO_REAL ( ITEM : D_NUMBER ) return REAL is
    NUMBER : REAL := 0.0 ;
    POWER : REAL := REAL ( BASE ) ;
  begin
    for I in reverse SIZE loop
      NUMBER := NUMBER + REAL ( ITEM.D_DIGITS( I )) / POWER ;
      POWER := POWER * REAL ( BASE ) ;
    end loop ;
    NUMBER := NUMBER * ( REAL( BASE ) ** ITEM.EXP) ;
    if ITEM.SIGN = MINUS then
      NUMBER := - NUMBER ;
    end if ;
    return NUMBER ;
  end D_NUMBER_TO_REAL ;

  function REAL_TO_D_NUMBER ( ITEM : REAL ) return D_NUMBER is
    NUMBER : D_NUMBER ;
    VALUE : REAL := abs ITEM ;
    REAL_BASE : REAL := REAL ( BASE ) ;
  begin
    if ITEM = 0.0 then
      return NUMBER ;
    end if ;
    if ITEM < 0.0 then
      NUMBER.SIGN := MINUS ;
    end if ;
    if VALUE >= 1.0 then
      while VALUE >= 1.0 loop
        NUMBER.EXP := NUMBER.EXP + 1 ;
        VALUE := VALUE / REAL_BASE ;
      end loop ;
    elsif VALUE < 1.0 / REAL_BASE then
      while VALUE < 1.0 / REAL_BASE loop
        NUMBER.EXP := NUMBER.EXP - 1 ;
        VALUE := VALUE * REAL_BASE ;
      end loop ;
    end if ;
    for I in reverse SIZE loop
      if VALUE >= 1.0 then
        VALUE := 1.0 - REAL'EPSILON ;
      end if ;
      NUMBER.D_DIGITS ( I ) := INTEGER ( VALUE * REAL_BASE - 0.5 ) ;
      if NUMBER.D_DIGITS ( I ) >= BASE then
        NUMBER.D_DIGITS ( I ) := BASE - 1 ;
      elsif NUMBER.D_DIGITS ( I ) < 0 then
        NUMBER.D_DIGITS ( I ) := 0 ;
      end if ;
      VALUE := VALUE * REAL_BASE - REAL ( NUMBER.D_DIGITS( I )) ;
      exit when VALUE <= REAL'EPSILON ;
      exit when SIZE'LAST - I >= REAL'DIGITS ;
    end loop ;
    return NORMALIZE ( NUMBER ) ;
  end REAL_TO_D_NUMBER ;


  procedure PUT ( FILE : in FILE_TYPE ;
                  ITEM : in D_NUMBER ;
                  FORE : in FIELD := DEFAULT_FORE ;
                  AFT  : in FIELD := DEFAULT_AFT ;
                  EXP  : in FIELD := DEFAULT_EXP ) is
  begin
    if FORE > 2 then
      TEXT_IO.PUT ( FILE , STRING'( 1..FORE-2 => ' ' ) ) ;
    end if ;
    if ITEM.SIGN = MINUS then
      TEXT_IO.PUT ( FILE , "-0." ) ;
    else
      TEXT_IO.PUT ( FILE , " 0." ) ;
    end if ;
    for I in reverse INTEGER ( SIZE'LAST ) - AFT .. INTEGER ( SIZE'LAST ) loop
      TEXT_IO.PUT ( FILE , INTEGER'IMAGE( ITEM.D_DIGITS(SIZE(I))) ( 2 .. 2 )) ;
    end loop ;
    TEXT_IO.PUT ( FILE , "E" ) ;
    INT_IO.PUT ( FILE , ITEM.EXP , WIDTH => 0 ) ;
  end PUT ;

  procedure PUT ( ITEM : in D_NUMBER ;
                  FORE : in FIELD := DEFAULT_FORE ;
                  AFT  : in FIELD := DEFAULT_AFT ;
                  EXP  : in FIELD := DEFAULT_EXP ) is
  begin
    if ITEM.SIGN = MINUS then
      TEXT_IO.PUT ( "-0." ) ;
    else
      TEXT_IO.PUT ( " 0." ) ;
    end if ;
    for I in reverse SIZE loop
      TEXT_IO.PUT ( INTEGER'IMAGE( ITEM.D_DIGITS( I )) ( 2 .. 2 )) ;
    end loop ;
    TEXT_IO.PUT ( "E" ) ;
    TEXT_IO.PUT ( INTEGER'IMAGE( ITEM.EXP )) ;
  end PUT ;

  procedure PUT ( TO   : out STRING ;
                  ITEM : in D_NUMBER ;
                  AFT  : in FIELD := DEFAULT_AFT ;
                  EXP  : in FIELD := DEFAULT_EXP ) is
    POS : POSITIVE := 4 ;
    LAST : POSITIVE ;
  begin
    if ITEM.SIGN = MINUS then
      TO ( 1..3 ) := ( "-0." ) ;
    else
      TO ( 1..3 ) := ( " 0." ) ;
    end if ;
    for I in reverse SIZE loop
      TO ( POS..POS ) := INTEGER'IMAGE( ITEM.D_DIGITS( I )) ( 2 .. 2 ) ;
      POS := POS + 1 ;
    end loop ;
    TO ( POS ) := 'E' ;
    LAST := INTEGER'IMAGE( ITEM.EXP )'LENGTH ;
    TO ( POS+1..POS+LAST ) := INTEGER'IMAGE( ITEM.EXP ) ;
    for I in POS+LAST+1..TO'LAST loop
      TO ( I ) := ' ' ;
    end loop ;
  end PUT ;

  procedure GET ( FILE : in FILE_TYPE ;
                  ITEM : out D_NUMBER ;
                  WIDTH : in FIELD := 0 ) is
    FROM : STRING ( 1..INTEGER(SIZE'LAST+15) ) ;
    LAST : POSITIVE ;
  begin
    TEXT_IO.GET_LINE ( FILE , FROM , LAST ) ;
    GET ( FROM(1..LAST) , ITEM , LAST ) ;
  end GET ;

  procedure GET ( ITEM : out D_NUMBER ;
                  WIDTH : in FIELD := 0 ) is
    FROM : STRING ( 1..INTEGER(SIZE'LAST+15) ) ;
    LAST : POSITIVE ;
  begin
    TEXT_IO.GET_LINE ( FROM , LAST ) ;
    GET ( FROM(1..LAST) , ITEM , LAST ) ;
  end GET ;

  procedure GET ( FROM : in STRING ;
                  ITEM : out D_NUMBER ;
                  LAST : out POSITIVE ) is
    STRING_POSITION : POSITIVE := 1 ;
    POINT_FOUND : BOOLEAN := FALSE ;
    AFTER_POINT : INTEGER := 0 ;
    DIGIT_POSITION : INTEGER := INTEGER(SIZE'LAST)+1 ;
    HOLD : D_NUMBER := CONSTANT_ZERO ;
    EXPONENT : INTEGER := 0 ;
    FROM_CHAR : CHARACTER := '#';
    EXPONENT_NEGATIVE : BOOLEAN := FALSE ;
  begin
    ITEM := CONSTANT_ZERO ;
-- skip leading blank and zero
    while STRING_POSITION <= FROM'LAST loop
      if FROM(STRING_POSITION) /= ' ' and
         FROM(STRING_POSITION) /= '0' then
        exit ;
      end if ;
      STRING_POSITION := STRING_POSITION + 1 ;
    end loop ;
-- pick up mantissa
    while STRING_POSITION <= FROM'LAST loop
      if FROM(STRING_POSITION) = '_' then
        null ;
      elsif FROM(STRING_POSITION) = '+' then
        null ;
      elsif FROM(STRING_POSITION) = '.' then
        POINT_FOUND := TRUE ;
      elsif FROM(STRING_POSITION) = '-' then
        HOLD.SIGN := MINUS ;
      elsif FROM(STRING_POSITION) = 'E' or
            FROM(STRING_POSITION) = 'e' then
         STRING_POSITION := STRING_POSITION + 1 ;
         exit ;
      elsif FROM(STRING_POSITION) = '0' then
        if DIGIT_POSITION /= INTEGER(SIZE'LAST)+1 and
           DIGIT_POSITION > 0 then
          DIGIT_POSITION := DIGIT_POSITION - 1 ;
          HOLD.D_DIGITS(SIZE(DIGIT_POSITION)) := 0 ;
        end if ;
        if POINT_FOUND then
          AFTER_POINT := AFTER_POINT + 1 ;
        end if ;
      elsif FROM(STRING_POSITION) = '1' then
        if DIGIT_POSITION > 0 then
          DIGIT_POSITION := DIGIT_POSITION - 1 ;
          HOLD.D_DIGITS(SIZE(DIGIT_POSITION)) := 1 ;
          if POINT_FOUND then
            AFTER_POINT := AFTER_POINT + 1 ;
          end if ;
        end if ;
      elsif FROM(STRING_POSITION) = '2' then
        if DIGIT_POSITION > 0 then
          DIGIT_POSITION := DIGIT_POSITION - 1 ;
          HOLD.D_DIGITS(SIZE(DIGIT_POSITION)) := 2 ;
          if POINT_FOUND then
            AFTER_POINT := AFTER_POINT + 1 ;
          end if ;
        end if ;
      elsif FROM(STRING_POSITION) = '3' then
        if DIGIT_POSITION > 0 then
          DIGIT_POSITION := DIGIT_POSITION - 1 ;
          HOLD.D_DIGITS(SIZE(DIGIT_POSITION)) := 3 ;
          if POINT_FOUND then
            AFTER_POINT := AFTER_POINT + 1 ;
          end if ;
        end if ;
      elsif FROM(STRING_POSITION) = '4' then
        if DIGIT_POSITION > 0 then
          DIGIT_POSITION := DIGIT_POSITION - 1 ;
          HOLD.D_DIGITS(SIZE(DIGIT_POSITION)) := 4 ;
          if POINT_FOUND then
            AFTER_POINT := AFTER_POINT + 1 ;
          end if ;
        end if ;
      elsif FROM(STRING_POSITION) = '5' then
        if DIGIT_POSITION > 0 then
          DIGIT_POSITION := DIGIT_POSITION - 1 ;
          HOLD.D_DIGITS(SIZE(DIGIT_POSITION)) := 5 ;
          if POINT_FOUND then
            AFTER_POINT := AFTER_POINT + 1 ;
          end if ;
        end if ;
      elsif FROM(STRING_POSITION) = '6' then
        if DIGIT_POSITION > 0 then
          DIGIT_POSITION := DIGIT_POSITION - 1 ;
          HOLD.D_DIGITS(SIZE(DIGIT_POSITION)) := 6 ;
          if POINT_FOUND then
            AFTER_POINT := AFTER_POINT + 1 ;
          end if ;
        end if ;
      elsif FROM(STRING_POSITION) = '7' then
        if DIGIT_POSITION > 0 then
          DIGIT_POSITION := DIGIT_POSITION - 1 ;
          HOLD.D_DIGITS(SIZE(DIGIT_POSITION)) := 7 ;
          if POINT_FOUND then
            AFTER_POINT := AFTER_POINT + 1 ;
          end if ;
        end if ;
      elsif FROM(STRING_POSITION) = '8' then
        if DIGIT_POSITION > 0 then
          DIGIT_POSITION := DIGIT_POSITION - 1 ;
          HOLD.D_DIGITS(SIZE(DIGIT_POSITION)) := 8 ;
          if POINT_FOUND then
            AFTER_POINT := AFTER_POINT + 1 ;
          end if ;
        end if ;
      elsif FROM(STRING_POSITION) = '9' then
        if DIGIT_POSITION > 0 then
          DIGIT_POSITION := DIGIT_POSITION - 1 ;
          HOLD.D_DIGITS(SIZE(DIGIT_POSITION)) := 9 ;
          if POINT_FOUND then
            AFTER_POINT := AFTER_POINT + 1 ;
          end if ;
        end if ;
      else
        TEXT_IO.PUT(FROM(STRING_POSITION));
        TEXT_IO.PUT_LINE ( " = unexpected character in mantissa " ) ;
      end if;
      STRING_POSITION := STRING_POSITION + 1 ;
    end loop ;
-- get exponent
    while STRING_POSITION <= FROM'LAST loop
      if FROM(STRING_POSITION) = '-' then
        EXPONENT_NEGATIVE := TRUE ;
      elsif FROM(STRING_POSITION) = '+' then
        null ;
      elsif FROM(STRING_POSITION) = '_' then
        null ;
      elsif FROM(STRING_POSITION) = ' ' then
        exit ;
      elsif FROM(STRING_POSITION) = '0' then
        EXPONENT := EXPONENT * 10 + 0 ;
      elsif FROM(STRING_POSITION) = '1' then
        EXPONENT := EXPONENT * 10 + 1 ;
      elsif FROM(STRING_POSITION) = '2' then
        EXPONENT := EXPONENT * 10 + 2 ;
      elsif FROM(STRING_POSITION) = '3' then
        EXPONENT := EXPONENT * 10 + 3 ;
      elsif FROM(STRING_POSITION) = '4' then
        EXPONENT := EXPONENT * 10 + 4 ;
      elsif FROM(STRING_POSITION) = '5' then
        EXPONENT := EXPONENT * 10 + 5 ;
      elsif FROM(STRING_POSITION) = '6' then
        EXPONENT := EXPONENT * 10 + 6 ;
      elsif FROM(STRING_POSITION) = '7' then
        EXPONENT := EXPONENT * 10 + 7 ;
      elsif FROM(STRING_POSITION) = '8' then
        EXPONENT := EXPONENT * 10 + 8 ;
      elsif FROM(STRING_POSITION) = '9' then
        EXPONENT := EXPONENT * 10 + 9 ;
      else
        TEXT_IO.PUT(FROM(STRING_POSITION));
        TEXT_IO.PUT_LINE ( " = unexpected character in exponent " ) ;
      end if ;
      STRING_POSITION := STRING_POSITION + 1 ;
    end loop ;
-- check result
    if EXPONENT_NEGATIVE then
      EXPONENT := - EXPONENT ;
    end if ;
    EXPONENT := EXPONENT + (INTEGER(SIZE'LAST)+1-DIGIT_POSITION) - AFTER_POINT;
    HOLD.EXP := EXPONENT ;
    if HOLD.D_DIGITS(SIZE'LAST) = 0 then
      ITEM := NORMALIZE ( HOLD ) ;
    else
      ITEM := HOLD ;
    end if ;
    LAST := STRING_POSITION ;
  end GET ;


  procedure DE_FLOAT ( X : D_NUMBER ;
                       SIGN : out INTEGER ;
                       MANTISSA : out D_NUMBER ;
                       EXPONENT : out INTEGER ) is
  begin
    MANTISSA := X ;
    MANTISSA.SIGN := PLUS ;
    MANTISSA.EXP := 0 ;
    if X.SIGN = PLUS then
      SIGN := 1 ;
    else
      SIGN := -1 ;
    end if ;
    EXPONENT := X.EXP ;
  end DE_FLOAT ;


  procedure RE_FLOAT ( X : out D_NUMBER ;
                       SIGN : INTEGER ;
                       MANTISSA : D_NUMBER ;
                       EXPONENT : INTEGER ) is
  begin
    X := MANTISSA ;
    if SIGN = 1 then
      X.SIGN := PLUS ;
    else
      X.SIGN := MINUS ;
    end if ;
    X.EXP := EXPONENT ;
  end RE_FLOAT ;

  function "rem" ( LEFT , RIGHT : D_NUMBER ) return D_NUMBER is
                 -- LEFT/RIGHT gives remainder less than RIGHT
    DENOMINATOR : D_NUMBER := abs LEFT ;
    DIVISOR : D_NUMBER := abs RIGHT ;
    REDUCER : D_NUMBER ;
    DIVISOR_MANTISSA : D_NUMBER ;
    EXPONENT : INTEGER ;
    MANTISSA : D_NUMBER ;
    SIGN : INTEGER ;

  begin
    if RIGHT.D_DIGITS(SIZE'LAST) = 0 then
      raise NUMERIC_ERROR ;
    elsif LEFT.D_DIGITS(SIZE'LAST) = 0 then
      return LEFT ;
    elsif LEFT = RIGHT then
      return CONSTANT_ZERO ;
    end if ;
    DE_FLOAT ( DIVISOR , SIGN , DIVISOR_MANTISSA , EXPONENT ) ;
    while DENOMINATOR >= DIVISOR loop
      DE_FLOAT ( DENOMINATOR , SIGN , MANTISSA , EXPONENT ) ;
      RE_FLOAT ( REDUCER , SIGN , DIVISOR_MANTISSA , EXPONENT ) ;
      if REDUCER > DENOMINATOR then
--        REDUCER := REDUCER * 0.1 ; -- needs to be exact, may underflow
        RE_FLOAT ( REDUCER , SIGN , DIVISOR_MANTISSA , EXPONENT - 1 ) ;
        if REDUCER.D_DIGITS(SIZE'LAST) = 0 then
          return CONSTANT_ZERO ;
        end if ;
      end if ;
      DENOMINATOR := DENOMINATOR - REDUCER ;
    end loop ;
    if LEFT.SIGN = MINUS then
      return - DENOMINATOR ; --the reduced denominator has become the remainder
    else
      return DENOMINATOR ;
    end if ;
  end "rem" ;

  function SQRT ( X : D_NUMBER ) return D_NUMBER is
    Y : D_NUMBER := X ;
    COUNT : INTEGER := 1 ;
  begin
    if X.SIGN = MINUS then
      raise NUMERIC_ERROR ;
    elsif X.D_DIGITS ( SIZE'LAST ) = 0 then
      return Y ;
    end if ;
    Y.EXP := X.EXP / 2 ;
    if abs X.EXP mod 2 = 1 then
      Y := Y * REAL_TO_D_NUMBER ( 3.162 ) ;
    end if ;
    Y := CONSTANT_HALF * ( Y + X / Y ) ;
    Y := CONSTANT_HALF * ( Y + X / Y ) ;
    while COUNT < INTEGER ( SIZE'LAST ) loop
      Y := CONSTANT_HALF * ( Y + X / Y ) ;
      COUNT := COUNT + COUNT ;
    end loop ;
    return Y ;
  end SQRT ;


  function LOG ( X : D_NUMBER ) return D_NUMBER is
--  X = MANTISSA * 10**EXPONENT
--  LN(X) := LOG10(X) * LN(10)
--           LOG10(X) = LOG10(10**EXPONENT)/LN(10) + LN(MANTISSA)
--
--  LN(F/2) = LN(F)-LN(2) reduction
--
    RESULT : D_NUMBER ;
    REDUCED_X : D_NUMBER ;
    XX : D_NUMBER ;
    SIGN : INTEGER ;
    MANTISSA : D_NUMBER ;
    EXPONENT : INTEGER ;
    I : INTEGER := ( ONE_OVER'LAST / 2 ) * 2 - 1 ;
  begin
    if X <= CONSTANT_ZERO then
      raise NUMERIC_ERROR ;
    end if ;
    DE_FLOAT ( X , SIGN , MANTISSA , EXPONENT ) ;
    if MANTISSA < CONSTANT_EIGHTH then
      REDUCED_X := MANTISSA * CONSTANT_EIGHT ;
      REDUCED_X := ( REDUCED_X - CONSTANT_ONE ) / ( REDUCED_X + CONSTANT_ONE ) ;
    elsif MANTISSA < CONSTANT_FOURTH then
      REDUCED_X := MANTISSA * CONSTANT_FOUR ;
      REDUCED_X := ( REDUCED_X - CONSTANT_ONE ) / ( REDUCED_X + CONSTANT_ONE ) ;
    elsif MANTISSA < CONSTANT_HALF then
      REDUCED_X := MANTISSA * CONSTANT_TWO ;
      REDUCED_X := ( REDUCED_X - CONSTANT_ONE ) / ( REDUCED_X + CONSTANT_ONE ) ;
    else
      REDUCED_X := ( MANTISSA - CONSTANT_ONE ) / ( MANTISSA + CONSTANT_ONE ) ;
    end if ;
    XX := REDUCED_X * REDUCED_X ;
    while I > 1 loop
      RESULT := ( ONE_OVER ( I ) + RESULT ) * XX ;
      I := I - 2 ;
    end loop ;
    RESULT := CONSTANT_TWO * ( CONSTANT_ONE + RESULT ) * REDUCED_X ;
    if MANTISSA < CONSTANT_EIGHTH then
      RESULT := RESULT - LOG_BASE_E_OF_2 - LOG_BASE_E_OF_2  - LOG_BASE_E_OF_2 ;
    elsif MANTISSA < CONSTANT_FOURTH then
      RESULT := RESULT - LOG_BASE_E_OF_2 - LOG_BASE_E_OF_2 ;
    elsif MANTISSA < CONSTANT_HALF then
      RESULT := RESULT - LOG_BASE_E_OF_2 ;
    end if ;
    RESULT := RESULT + INTEGER_TO_D_NUMBER ( EXPONENT ) * LOG_BASE_E_OF_10 ;
    return RESULT ;
  end LOG ;

  function EXP ( X : D_NUMBER ) return D_NUMBER is
    RESULT : D_NUMBER := CONSTANT_ZERO ;
    REDUCED_X : D_NUMBER := abs X ;
    INTEGER_PART : INTEGER ;
    POWER : D_NUMBER := CONSTANT_ONE ;
  begin
    if abs X >= CONSTANT_ONE then
      INTEGER_PART := D_NUMBER_TO_INTEGER ( abs X ) ;
      POWER := CONSTANT_NATURAL_E ** INTEGER_PART ;
      REDUCED_X := ( abs X ) rem CONSTANT_ONE ; -- get fractional part
    end if ;
    REDUCED_X := CONSTANT_HALF * REDUCED_X ;
    for I in reverse 1 .. FACTORIAL'LAST loop
      RESULT := ( FACTORIAL ( I ) + RESULT ) * REDUCED_X ;
    end loop ;
    RESULT := RESULT + CONSTANT_ONE ;
    RESULT := POWER * RESULT * RESULT ;
    if X.SIGN = MINUS then
      return CONSTANT_ONE / RESULT ;
    else
      return RESULT ;
    end if ;
  end EXP ;

  function ATAN2 ( Y , X : D_NUMBER ) return D_NUMBER is
    Z : D_NUMBER := CONSTANT_ZERO ;
    A : D_NUMBER ;
    A_NORM : D_NUMBER ;
    A_2 : D_NUMBER ;
    REDUCE : D_NUMBER := CONSTANT_ONE ;
    F_REDUCE_CASE : INTEGER := 0 ;
    F_CASE : array ( 0..3 ) of D_NUMBER := ( CONSTANT_ZERO ,
        CONSTANT_PI_OVER_SIX , CONSTANT_PI_OVER_TWO , CONSTANT_PI_OVER_THREE );
  begin
    if X.D_DIGITS(SIZE'LAST) = 0 and
       Y.D_DIGITS(SIZE'LAST) = 0 then
      raise NUMERIC_ERROR ;  -- by definition
    elsif abs Y > abs X then
      A := X / Y ;
      F_REDUCE_CASE := 2 ;
    else
      A := Y / X ;
    end if ;
    A_NORM := abs A ;
    if A_NORM > CONSTANT_ATAN_BREAK_1 then
      A_NORM := ( A_NORM * CONSTANT_SQRT_3 - CONSTANT_ONE ) /
                ( A_NORM + CONSTANT_SQRT_3 ) ;
      F_REDUCE_CASE := F_REDUCE_CASE + 1 ;
    end if ;
    if abs A_NORM > CONSTANT_ATAN_BREAK_2 then
      A_NORM := ( SQRT ( CONSTANT_ONE + A_NORM * A_NORM ) - CONSTANT_ONE ) /
                  A_NORM ;
      REDUCE := CONSTANT_TWO ;
    end if ;
    A_2 := A_NORM * A_NORM ;
    if abs A_2 < DIGITS_EPSILON then
      Z := A_2 ;
    else
      for I in reverse ATAN_C'RANGE loop
        Z := ( ATAN_C( I ) + Z) * A_2 ;
      end loop ;
      Z := ( CONSTANT_ONE + Z ) * A_NORM ;
    end if ;
    if F_REDUCE_CASE > 1 then
      Z := - Z ;
    end if ;
    Z := F_CASE ( F_REDUCE_CASE ) + REDUCE * Z ;
    if A.SIGN = MINUS then
      Z := - Z ;
    end if ;
    if abs Y > abs X then
      if Y.SIGN = MINUS then
        Z := - ( CONSTANT_PI_OVER_TWO + Z ) ;
      else
        Z := CONSTANT_PI_OVER_TWO - Z ;
      end if ;
    else
      if X.SIGN = MINUS and
         Y.SIGN = MINUS then
        Z := - CONSTANT_PI + Z ;
      elsif X.SIGN = MINUS and
         Y > CONSTANT_ZERO then
        Z := CONSTANT_PI + Z ;
      end if ;
    end if ;
    return Z ;
  end ATAN2 ;

  function SIN ( X : D_NUMBER ) return D_NUMBER is
    RESULT : D_NUMBER := CONSTANT_ZERO ;
    REDUCED_X : D_NUMBER ;
    XX : D_NUMBER ;
    NEGATE : BOOLEAN := TRUE ;
    I : INTEGER ;
  begin
    if abs X > CONSTANT_PI then
      REDUCED_X := X rem CONSTANT_TWO_PI ;
    else
      REDUCED_X := X ;
    end if ;
    if REDUCED_X > CONSTANT_PI then
      REDUCED_X := REDUCED_X - CONSTANT_TWO_PI ;
    elsif REDUCED_X < - CONSTANT_PI then
      REDUCED_X := REDUCED_X + CONSTANT_TWO_PI ;
    end if ;
    if REDUCED_X > CONSTANT_PI_OVER_TWO then
      REDUCED_X := CONSTANT_PI - REDUCED_X ;
    elsif REDUCED_X < - CONSTANT_PI_OVER_TWO then
      REDUCED_X := - ( CONSTANT_PI + REDUCED_X ) ;
    end if ;
    XX := REDUCED_X * REDUCED_X ;
    I := ( ( FACTORIAL'LAST / 4 ) * 4 ) - 1 ;
    while I > 1 loop
      if NEGATE then
        RESULT := ( - FACTORIAL ( I ) + RESULT ) * XX ;
      else
        RESULT := ( FACTORIAL ( I ) + RESULT ) * XX ;
      end if ;
      NEGATE := not NEGATE ;
      I := I - 2 ;
    end loop ;
    return ( CONSTANT_ONE + RESULT ) * REDUCED_X ;
  end SIN ;

  function NATURAL_E return D_NUMBER is
  begin
    return CONSTANT_NATURAL_E ; -- computed below after package begin
  end NATURAL_E ;

  function PI return D_NUMBER is
  begin
    return CONSTANT_PI ; -- computed below after package begin
  end PI ;


  function LN_10 return D_NUMBER is
  begin
    return LOG_BASE_E_OF_10 ; -- computed below after package begin
  end LN_10 ;

  function LN_2 return D_NUMBER is
  begin
    return LOG_BASE_E_OF_2 ; -- computed below after package begin
  end LN_2 ;

  function NORMALIZE ( ITEM : D_NUMBER ) return D_NUMBER is
    RESULT : D_NUMBER ;
    OFFSET : INTEGER ;
  begin
    RESULT.EXP := ITEM.EXP ;
    RESULT.SIGN := ITEM.SIGN ;
    if ITEM.D_DIGITS ( SIZE'LAST ) = 0 then
      OFFSET := 0 ;
      for I in reverse SIZE loop
        exit when ITEM.D_DIGITS ( I ) /= 0 ;
        OFFSET := OFFSET + 1 ;
        RESULT.EXP := RESULT.EXP - 1 ;
      end loop ;
      if OFFSET = INTEGER ( SIZE'LAST + 1 ) then  -- natural zero
        RESULT.SIGN := PLUS ;
        RESULT.EXP := 0 ;
      else
        for I in SIZE ( OFFSET ) .. SIZE'LAST loop
          RESULT.D_DIGITS ( I ) := ITEM.D_DIGITS ( I - SIZE( OFFSET )) ;
        end loop ;
      end if ;
    else  -- already normalized
      RESULT := ITEM ;
    end if ;
    return RESULT ;
  end NORMALIZE ;

begin

  if SIZE'FIRST /= 0 then
    raise PROGRAM_ERROR ;
  end if ;

-- basic digits constants

  CONSTANT_ZERO := DIGITS_ZERO ;
  CONSTANT_ONE := DIGITS_ONE ;
  CONSTANT_TWO := CONSTANT_ONE + CONSTANT_ONE ;
  CONSTANT_THREE := CONSTANT_TWO + CONSTANT_ONE ;
  CONSTANT_FOUR := CONSTANT_THREE + CONSTANT_ONE ;
  CONSTANT_SIX := CONSTANT_THREE + CONSTANT_THREE ;
  CONSTANT_EIGHT := CONSTANT_FOUR + CONSTANT_FOUR ;
  CONSTANT_TEN := CONSTANT_EIGHT + CONSTANT_TWO ;
  CONSTANT_TWENTY_FOUR := CONSTANT_EIGHT * CONSTANT_THREE ;
  CONSTANT_HALF := CONSTANT_ONE / CONSTANT_TWO ;
  CONSTANT_FOURTH := CONSTANT_ONE / CONSTANT_FOUR ;
  CONSTANT_EIGHTH := CONSTANT_ONE / CONSTANT_EIGHT ;
  CONSTANT_EPSILON.EXP := - ( INTEGER ( SIZE'LAST ) - 1 ) ;
  CONSTANT_EPSILON.D_DIGITS ( SIZE'LAST ) := 1 ;
  CONSTANT_SQRT_3 := REAL_TO_D_NUMBER ( 1.732050807568877293527446341505872) ;
  while abs ( CONSTANT_SQRT_3 * CONSTANT_SQRT_3 - CONSTANT_THREE ) >
                                  CONSTANT_FOUR * CONSTANT_EPSILON loop
    CONSTANT_SQRT_3 := CONSTANT_HALF * ( CONSTANT_SQRT_3 + CONSTANT_THREE /
                                         CONSTANT_SQRT_3 ) ;
  end loop ;
  CONSTANT_ONE_OVER_SQRT_3 := CONSTANT_ONE / CONSTANT_SQRT_3 ;
  CONSTANT_ATAN_BREAK_1 := CONSTANT_TWO - CONSTANT_SQRT_3 ;
  SQRT_ATAN_BREAK_1 := REAL_TO_D_NUMBER (
                                      0.5176380902050415246977976752480966 ) ;
  while abs ( SQRT_ATAN_BREAK_1 * SQRT_ATAN_BREAK_1 - CONSTANT_ATAN_BREAK_1 ) >
                                        CONSTANT_FOUR * CONSTANT_EPSILON loop
    SQRT_ATAN_BREAK_1 := CONSTANT_HALF * ( SQRT_ATAN_BREAK_1 +
                         CONSTANT_ATAN_BREAK_1 / SQRT_ATAN_BREAK_1 ) ;
  end loop ;
  CONSTANT_ATAN_BREAK_2 := ( CONSTANT_TWO * SQRT_ATAN_BREAK_1 - CONSTANT_ONE )
                              / CONSTANT_ATAN_BREAK_1 ;

-- compute reciprocal factorials for use by functions

  FACTORIAL ( 0 ) := CONSTANT_ONE ;
  FACTORIAL ( 1 ) := CONSTANT_ONE ;
  for I in 2..FACTORIAL'LAST loop
    FACTORIAL ( I ) := FACTORIAL ( I - 1 ) / INTEGER_TO_D_NUMBER ( I ) ;
  end loop ;

-- compute reciprocals for use by LOG
  for I in 1..ONE_OVER'LAST loop
    ONE_OVER ( I ) := CONSTANT_ONE / INTEGER_TO_D_NUMBER ( I ) ;
  end loop ;

-- compute ATAN coefficients
  declare
    ATAN_COUNT : INTEGER := - 3 ;
  begin
    for I in ATAN_C'RANGE loop
      ATAN_C(I) := CONSTANT_ONE / INTEGER_TO_D_NUMBER ( ATAN_COUNT ) ;
      if ATAN_COUNT < 0 then
        ATAN_COUNT := - ATAN_COUNT + 2 ;
      else
        ATAN_COUNT := - ( ATAN_COUNT + 2 ) ;
      end if ;
    end loop ;
  end ;

  declare -- compute PI   6 * ARCTAN ( 1/SQRT(3) ) -- SIZE'LAST*4+3 terms
--                       12 * ARCTAN ( 2-SQRT(3) ) -- SIZE'LAST*2+  terms
--                       24 * ARCTAN ( (2*SQRT(2-SQRT(3))-1)/(2-SQRT(3) )
    Y : D_NUMBER := CONSTANT_ZERO ;
    A : D_NUMBER := CONSTANT_ATAN_BREAK_2 ;
    A_2 : D_NUMBER ;
  begin
    A_2 := A * A ;
    for J in reverse ATAN_C'RANGE loop
      Y := ( ATAN_C ( J ) + Y ) * A_2 ;
    end loop ;
    CONSTANT_PI := CONSTANT_TWENTY_FOUR * ( CONSTANT_ONE + Y ) * A ;
  end ;

  declare -- compute NATURAL_E
    RESULT : D_NUMBER := CONSTANT_ZERO ;
  begin
    for I in reverse SIZE loop
      RESULT := RESULT + FACTORIAL ( INTEGER ( I ) ) ;
    end loop ;
    CONSTANT_NATURAL_E := RESULT ;
  end ;

-- compute constants that depend on functions above
  CONSTANT_TWO_PI := CONSTANT_TWO * CONSTANT_PI ;
  CONSTANT_PI_OVER_TWO := CONSTANT_HALF * CONSTANT_PI ;
  CONSTANT_PI_OVER_THREE := CONSTANT_PI / CONSTANT_THREE ;
  CONSTANT_PI_OVER_FOUR := CONSTANT_PI / CONSTANT_FOUR ;
  CONSTANT_PI_OVER_SIX := CONSTANT_PI / CONSTANT_SIX ;

  declare -- get  LN(2) and LN(10)
    RESULT : D_NUMBER := CONSTANT_ZERO ;
    X_2 : D_NUMBER := CONSTANT_TWO /  CONSTANT_NATURAL_E ;
    X_10 : D_NUMBER := CONSTANT_TEN / ( CONSTANT_NATURAL_E *
                           CONSTANT_NATURAL_E  * CONSTANT_NATURAL_E ) ;
    REDUCED_X : D_NUMBER ;
    XX : D_NUMBER ;
    I : INTEGER ;
  begin
    REDUCED_X := ( X_2 - CONSTANT_ONE ) / ( X_2 + CONSTANT_ONE ) ;
    XX := REDUCED_X * REDUCED_X ;
    I := ( ONE_OVER'LAST / 2 ) * 2 - 1 ;
    while I > 1 loop
      RESULT := ( ONE_OVER ( I ) + RESULT ) * XX ;
      I := I - 2 ;
    end loop ;
    RESULT := CONSTANT_TWO * ( CONSTANT_ONE + RESULT ) * REDUCED_X ;
    LOG_BASE_E_OF_2 := CONSTANT_ONE + RESULT ;
--
    RESULT := CONSTANT_ZERO ;
    X_10 := X_10 * CONSTANT_TWO ;
    REDUCED_X := ( X_10 - CONSTANT_ONE ) / ( X_10 + CONSTANT_ONE ) ;
    XX := REDUCED_X * REDUCED_X ;
    I := ( ONE_OVER'LAST / 4 ) * 2 - 1 ;
    while I > 1 loop
      RESULT := ( ONE_OVER ( I ) + RESULT ) * XX ;
      I := I - 2 ;
    end loop ;
    RESULT := CONSTANT_TWO * ( CONSTANT_ONE + RESULT ) * REDUCED_X ;
    LOG_BASE_E_OF_10 := CONSTANT_THREE + RESULT - LOG_BASE_E_OF_2 ;
  end ;

end GENERIC_DIGITS_ARITHMETIC ;