home *** CD-ROM | disk | FTP | other *** search
/ Programmer's ROM - The Computer Language Library / programmersrom.iso / ada / piwg / z000002.ada < prev    next >
Encoding:
Text File  |  1988-05-03  |  3.1 KB  |  119 lines
  1.  
  2. package REFUNCT is
  3.  
  4. -- stripped down real functions package
  5. -- this should work on any compiler on any computer 
  6.  
  7.   function EXP10 ( X : FLOAT ) return FLOAT ;
  8.  
  9.   function "**" ( X , Y : FLOAT ) return FLOAT ;
  10.  
  11.   function LOG10 ( X : FLOAT ) return FLOAT ;
  12.  
  13.   function SQRT ( X : FLOAT ) return FLOAT ;
  14.  
  15. end REFUNCT ;
  16.  
  17. package body REFUNCT is
  18.  
  19.   function EXP10 ( X : FLOAT ) return FLOAT is
  20.     C1 : constant FLOAT := 1.15129277603 ;
  21.     C2 : constant FLOAT := 0.66273088429 ;
  22.     C3 : constant FLOAT := 0.25439357484 ;
  23.     C4 : constant FLOAT := 0.07295173666 ;
  24.     C5 : constant FLOAT := 0.01742111988 ;
  25.     C6 : constant FLOAT := 0.00255491796 ;
  26.     C7 : constant FLOAT := 0.00093264267 ;
  27.     X1 : FLOAT ;
  28.     Y : FLOAT ;
  29.     TEN_PWR : FLOAT ;
  30.   begin
  31.     X1 := abs ( X ) ;
  32.     TEN_PWR := 1.0 ;
  33.     while X1 >= 1.0 loop
  34.       TEN_PWR := TEN_PWR * 10.0 ;
  35.       X1 := X1 - 1.0 ;
  36.     end loop ;
  37.     Y := 1.0 + ( C1 +( C2 +( C3 +( C4 +( C5 +( C6 + C7 * X1 ) * X1) * X1) * X1
  38.        ) * X1) * X1) * X1 ;
  39.     Y := Y * Y * TEN_PWR ;
  40.     if X < 0.0 then
  41.       Y := 1.0 / Y ;
  42.     end if ;
  43.     return Y ;
  44.   end EXP10 ;
  45.  
  46.   function "**" ( X , Y : FLOAT ) return FLOAT is
  47.   begin
  48.     if X = 0.0 then
  49.       return 0.0 ;
  50.     elsif Y = 0.0 then
  51.       return 1.0 ;
  52.     else
  53.       return EXP10 ( Y * LOG10( X )) ;
  54.     end if ;
  55.   end "**" ;
  56.  
  57.   function LOG10 ( X : FLOAT ) return FLOAT is
  58.     C1 : constant FLOAT := 0.868591718 ;
  59.     C3 : constant FLOAT := 0.289335524 ;
  60.     C5 : constant FLOAT := 0.177522071 ;
  61.     C7 : constant FLOAT := 0.094376476 ;
  62.     C9 : constant FLOAT := 0.191337714 ;
  63.     C_R10 : constant FLOAT := 3.1622777 ;
  64.     Y : FLOAT ;
  65.     X_NORM : FLOAT ;
  66.     X_LOG : FLOAT ;
  67.     FRAC : FLOAT ;
  68.     FRAC_2 : FLOAT ;
  69.   begin
  70.     X_LOG := 0.5 ;
  71.     X_NORM := X ;
  72.     if X >= 10.0 then
  73.       while X_NORM >= 10.0  -- REDUCE TO 1.0 .. 10.0
  74.           loop
  75.         X_LOG := X_LOG + 1.0 ;
  76.         X_NORM := X_NORM * 0.1 ;
  77.       end loop ;
  78.     else
  79.       while X_NORM < 1.0  -- REDUCE TO 1.0 .. 10.0
  80.           loop
  81.         X_LOG := X_LOG - 1.0 ;
  82.         X_NORM := X_NORM * 10.0 ;
  83.       end loop ;
  84.     end if ;
  85.     FRAC := ( X_NORM - C_R10 ) / ( X_NORM + C_R10 ) ;
  86.     FRAC_2 := FRAC * FRAC ;
  87.     Y := ( C1 +( C3 +( C5 +( C7 + C9 * FRAC_2 ) * FRAC_2) * FRAC_2) * FRAC_2)
  88.         * FRAC ;
  89.     return Y + X_LOG ;
  90.   end LOG10 ;
  91.  
  92.   function SQRT ( X : FLOAT ) return FLOAT is
  93.     Y , ROOT_PWR , X_NORM : FLOAT ;
  94.     A : constant FLOAT := 2.1902 ;
  95.     B : constant FLOAT := - 3.0339 ;
  96.     C : constant FLOAT := 1.5451 ;
  97.   begin
  98.     X_NORM := X ;
  99.     ROOT_PWR := 1.0 ;
  100.     if X > 1.0 then -- REDUCE TO 0.25 .. 1.0
  101.       while X_NORM > 1.0 loop
  102.         ROOT_PWR := ROOT_PWR * 2.0 ;
  103.         X_NORM := X_NORM * 0.25 ;
  104.       end loop ;
  105.     else
  106.       while X_NORM < 0.25 loop
  107.         ROOT_PWR := ROOT_PWR * 0.5 ;
  108.         X_NORM := X_NORM * 4.0 ;
  109.       end loop ;
  110.     end if ;
  111.     Y := A + B / ( C + X_NORM ) ;
  112.     Y := 0.5 * ( Y + X_NORM / Y ) ;
  113.     Y := 0.5 * ( Y + X_NORM / Y ) ;
  114.     Y := Y * ROOT_PWR ;
  115.     return Y ;
  116.   end SQRT ;
  117.  
  118. end REFUNCT ;
  119.