/* --------------------------------------------------------------------
--
--
-- Copyright©2002 by the Institute of Electrical and Electronics Engineers, Inc.
-- Three Park Avenue
-- New York, NY 10016-5997, USA
-- All rights reserved.
-- 
-- This document is an unapproved draft of a proposed IEEE Standard. As such,
-- this document is subject to change. USE AT YOUR OWN RISK! Because this
--  is an unapproved draft, this document must not be utilized for any 
-- conformance/compliance purposes. Permission is hereby granted for IEEE 
-- Standards Committee participants to reproduce this document for purposes 
-- of IEEE standardization activities only. Prior to submitting this document 
-- to another standards development organization for standardization 
-- activities, permission must first be obtained from the Manager, Standards 
-- Licensing and Contracts, IEEE Standards Activities Department. Other 
-- entities seeking permission to reproduce this document, in whole or in 
-- part, must obtain permission from the Manager, Standards Licensing and 
-- Contracts, IEEE Standard Activities Department.
--
-- IEEE Standards Activities Department
-- Standards Licensing and Contracts
-- 445 Hoes Lane, P.O. Box 1331
-- Piscataway, NJ 08855-1331, USA
--
--   Title     :  fphdl64_real_functions_base.inc < IEEE std # 1076.3 >
--
--   Developers:  VHDL Synthesis working group, PAR 1076.3
--
--   Purpose   :  Base functions for real based 64 bit floating point 
--
--   Limitation:
--
-- --------------------------------------------------------------------
--   Last Modified $Date: 2003-01-15 09:51:12-05 $  $Id: fphdl64_real_functions_base.inc,v 1.4 2003-01-15 09:51:12-05 bishop Exp $
--   $Log: fphdl64_real_functions_base.inc,v $
--   Revision 1.4  2003-01-15 09:51:12-05  bishop
--   Updated header
--
--   Revision 1.3  2003-01-15 09:22:43-05  bishop
--   1/15/2003 update
--
--   Revision 1.2  2003-01-03 13:02:15-05  bishop
--   Added trig functions
--
--   Revision 1.1  2003-01-02 09:41:11-05  bishop
--   Initial revision
--
-- --------------------------------------------------------------------
*/
// by Kurt Baty

//  -- purpose: Generates the necessary boundary number
//  -- Synthisable (64 bit only)
function [11:-52] fp_gen_number;
  input [3:0] fpstate; // -- type of number

  reg [11:-52] fpresult;

begin // -- function fp_gen_number
    fpresult = 0;
    case(fpstate)
	//  type valid_fpstate
	//	 nan = 0,           -- Signaling NaN
	//       quiet_nan = 1,     -- Quiet NaN
	//       neg_inf = 2,       -- Negative infinity
	//       neg = 3,           -- negative normalized nonzero
	//       neg_denormal = 4,  -- negative denormalized
	//       neg_zero = 5,      -- -0
	//       zero = 6,          -- +0
	//       denormal = 7,      -- Positive denormalized
	//       normal = 8,        -- positive normalized nonzero
	//       infinity = 9);     -- positive infinity
	//  -- valid fp states
    4'd0: begin // nan
        fpresult[11-1:0] = {11{1'b1}};
        // -- Exponent all "1"
        fpresult[-1] = 1'b1; // -- MSB of Fraction "1"
        // -- Note: From W. Khan "IEEE Standard 754 for Binary Floating Point"
        // -- The difference between a signaling NAN and a quiet NAN is that
        // -- the MSB of the Fraction is a "1" in a Signaling NAN, and is a
        // -- "0" in a quiet NAN.
    end
    4'd1: begin // quiet_nan
        fpresult[11-1:0] = {11{1'b1}};
        // -- Exponent all "1"
        fpresult[-52] = 1'b1; // -- LSB of Fraction "1" (Could have been any bit)
    end
    4'd2: begin // neg_inf
        fpresult[11-1:0] = {11{1'b1}};
        // -- Exponent all "1"
        fpresult[11] = 1'b1; // -- Sign bit to "1"
    end
    4'd9: begin // infinity
        fpresult[11-1:0] = {11{1'b1}};
        // -- Exponent all "1"
        // -- Sign bit is "0"
    end
    4'd8: begin // normal
        fpresult[11-2:0] = {11{1'b1}};
        // -- return 1.0
    end
    4'd3: begin // normal
        fpresult[11-2:0] = {11{1'b1}};
        fpresult[11] = 1'b1; // -- Sign bit to "1"
        // -- return -1.0
    end
    4'd6: begin // zero
        fpresult = 0;
    end
    4'd5: begin // neg_zero
        fpresult[11] = 1'b1; // -- Sign bit to "1"
    end
    default: begin // all other cases
        fpresult = 0;
    end
    endcase
    fp_gen_number = fpresult; // return fpresult;
end
endfunction // fp_gen_number


//  -- Returns the class which X falls into
//  -- Synthisable
function [3:0] fp_class;
        input [11:-52] x;   //floating point input
begin   // -- fp_class
        //  type valid_fpstate
        //       nan = 0,           -- Signaling NaN
        //       quiet_nan = 1,     -- Quiet NaN
        //       neg_inf = 2,       -- Negative infinity
        //       neg = 3,           -- negative normalized nonzero
        //       neg_denormal = 4,  -- negative denormalized
        //       neg_zero = 5,      -- -0
        //       zero = 6,          -- +0
        //       denormal = 7,      -- Positive denormalized
        //       normal = 8,        -- positive normalized nonzero
        //       infinity = 9);     -- positive infinity
        //  -- valid fp states
    // -- Special cases, check for illegal number
    if (&x[11-1:0]) begin
      // -- Exponent is all "1".
      if (|x[-1:-52]) begin
        // -- Fraction must be all "0" or this is not a number.
        if (x[-1])       // -- From "W. Khan - IEEE standard
                         // -- 754 binary FP"
            fp_class = 0;  // return nan; -- Signaling nan (Not a number)
        else
            fp_class = 1;  // return quiet_nan;
      end
      // -- Check for infinity
      else begin
        if (x[11] == 0)
          fp_class = 9;  // return infinity;   -- Positive infinity
        else
          fp_class = 2;  // return neg_inf;    -- Negative infinity
      end
    end  
    // -- check for "0"
    else begin
      if (~|x[11-1:0]) begin
        // -- Exponent is all "0"
        if (~|x[-1:-52]) begin
          // = '0' then                        -- Fraction is all "0"
          if (x[11] == 0)
            fp_class = 6;  // return zero;     -- Zero
          else
            fp_class = 5;  // return neg_zero;
        end
        else begin
          if (x[11] == 0)
            fp_class = 7;  // return denormal; -- Denormal number (ieee extended fp)
          else
            fp_class = 4;  // return neg_denormal;
        end
      end  
      else begin
        if (x[11] == 0)
          fp_class = 8;  // return normal;     -- Normal FP number
        else
          fp_class = 3;  // return neg;
      end
    end  
end
endfunction // fp_class


//  -- Arithmetic functions
//  -- Synthisable
function [63:0] absolute;
input [63:0] arg; // -- floating point input
begin
    absolute = {1'b0,arg[62:0]};
end
endfunction // absolute


//  -- IEEE 754 "negative" function
//  -- Synthisable
function [63:0] negative;
input [63:0] arg; // -- floating point input
begin
    negative = {~arg[63],arg[62:0]};
end
endfunction // negative

  
//  -- NOT Synthisable
function [63:0] addition;
   input [63:0] l, r;
   real real_l,real_r,real_res;
begin
   real_l = $bitstoreal(l);
   real_r = $bitstoreal(r);
   real_res = real_l + real_r;
   addition= $realtobits(real_res);
end
endfunction // addition

//  -- NOT Synthisable
function [63:0] subtract;
   input [63:0] l, r;
   real real_l,real_r,real_res;
begin
   real_l = $bitstoreal(l);
   real_r = $bitstoreal(r);
   real_res = real_l - real_r;
   subtract = $realtobits(real_res);
end
endfunction // subtract

//  -- NOT Synthisable
function [63:0] multiply;
   input [63:0] l, r;
   real real_l,real_r,real_res;
begin
   real_l = $bitstoreal(l);
   real_r = $bitstoreal(r);
   real_res = real_l * real_r;
   multiply = $realtobits(real_res);
end
endfunction // multiply

//  -- NOT Synthisable
function [63:0] divide;
   input [63:0] l, r;
   real real_l,real_r,real_res;
begin
   real_l = $bitstoreal(l);
   real_r = $bitstoreal(r);
   real_res = real_l / real_r;
   divide = $realtobits(real_res);
end
endfunction // divide

//  -- NOT Synthisable
function [63:0] remainder;
   input [63:0] l, r;
   real real_l,real_r,real_res;
   reg signed [63:0] int_part;
begin
   real_l = $bitstoreal(l);
   real_r = $bitstoreal(r);
   real_res = real_l / real_r;
   int_part = real_res;
   real_res = real_l - real_r * int_part;
   remainder = $realtobits(real_res);
end
endfunction // remainder

//  -- NOT Synthisable
function [63:0] square_root;
   input [63:0] arg;
   real	real_arg,root,root_last;
begin
  if(&(arg | ~arg) === 1'b1) begin
    real_arg = $bitstoreal(absolute(arg));
    root_last = 1;
    root = (1 + real_arg) / 2;
    while(root_last != root) begin
        root_last = root;
	root = (root_last + (real_arg/root_last))/2;
    end	
    square_root = $realtobits(root); 
  end
  else
    square_root = 64'bx;
end
endfunction // square_root


//  -- NOT Synthisable
function [63:0] logarithm;
   input signed [31:0] base;
   input [63:0] arg;
   real real_base,real_log;
   integer i;
begin
  if(&(arg | ~arg) === 1'b1) begin
      i = base;
      real_base = $itor(i);
      real_log = $bitstoreal(natural_log(arg))
               / $bitstoreal(natural_log($realtobits(real_base)));    
      logarithm = $realtobits(real_log);    
  end
  else
    logarithm = 64'bx;
end
endfunction // logarithm

 
//  -- NOT Synthisable
function [63:0] natural_log;
   input [63:0] arg;
   reg [3:0] arg_fptype;
   real real_arg,term,term1,ln,ln_last;
   integer i;
   reg signed [12-1:0] signed_exp;
   reg [52-1:0] sargfract;
   reg [63:0] ln_bits;
begin
  if(&(arg | ~arg) === 1'b1) begin
    arg_fptype = fp_class(arg);
    case(arg_fptype)
    4'd2,4'd3,4'd4: begin
	natural_log = fp_gen_number(4'd1); // ln of a negative is NAN
    end
    4'd5,4'd6: begin
	natural_log = fp_gen_number(4'd2); // ln of a zero is neg_inf
    end
    4'd7,4'd8: begin
        if(arg_fptype == 4'd7) begin
            signed_exp = {~arg[63-1],~arg[63-1],arg[63-2:63-12+1]} + 2'd2;
            for(i=1;i<=52;i=i+1) begin
                if(arg[52-i]) begin
                    signed_exp = signed_exp - i;
                    sargfract = arg[52-1:0] << i;
                    i = 53;
                end
            end
            real_arg = $bitstoreal({12'h3ff,sargfract});
        end
        else begin
            signed_exp = {~arg[63-1],~arg[63-1],arg[63-2:63-12+1]} + 1'b1;
            real_arg = $bitstoreal({12'h3ff,arg[52-1:0]});
        end
        ln_last = 0;
        term = (real_arg - 1.0) / (real_arg + 1);
        ln = term;
        term1 = term * term;
        i=1;
        while(ln_last != ln) begin
            ln_last = ln;
            term = i * term * term1 / (i+2);
            ln = ln_last + term;
            i = i+2;
        end
        ln = ln * 2.0 + signed_exp * $bitstoreal(64'h3fe6_2e42_fefa_39ef);
        natural_log = $realtobits(ln);
    end
    4'd9: begin
	natural_log = fp_gen_number(4'd9); // ln of a infinity is infinity
    end
    default: begin
	natural_log = fp_gen_number(4'd1); // ln is NAN
    end
    endcase
  end
  else
    natural_log = 64'bx;
end
endfunction // natural_log


//  -- NOT Synthisable
function [63:0] power_of_fp;
   input [63:0] l, r;
   reg [3:0] l_fptype,r_fptype;
   reg signed [31:0] signed_reg_r;
   integer i,integer_r;
   reg	r_is_an_integer;
   real real_l,real_r,log2_abs_l,log2_result;
   reg signed [31:0] signed_reg_log2_result;
   integer integer_log2_result;
   real log2_result_fraction,term,term1,pow,pow_last;
   reg signed [11-1:0] signed_exp;
   reg [52:-3] fraction;
   reg [63:0] pow_bits;
   reg [3:0] pow_fptype;
begin
  if(    (&(l | ~l) === 1'b1)
      && (&(r | ~r) === 1'b1)) begin
    l_fptype = fp_class(l);
    r_fptype = fp_class(r);
    signed_reg_r = fp64_to_integer_func(r);
    integer_r = signed_reg_r;
    real_l = $bitstoreal(l);
    real_r = $bitstoreal(r);
    r_is_an_integer = (real_r == integer_r);
    casex({r_is_an_integer,l_fptype,r_fptype})
    9'hx0x,9'hxx0,
    9'hx1x,9'hxx1,        // l**NAN, NAN**r is NAN
    9'hx37,9'hx47,
    9'hx2x,
    9'h03x,9'h04x: begin  // a negative ** non-integer is NAN
      power_of_fp = fp_gen_number(4'd1); // return NAN
    end
    9'hx52,9'hx62,        // 0**neg_infinity is 0.0
    9'hx53,9'hx63,
    9'hx54,9'hx64,        // 0**-real is 0.0
    9'hx55,9'hx65,
    9'hx56,9'hx66,        // 0**0 is 0.0
    9'hx57,9'hx67,
    9'hx58,9'hx68,        // 0**+real is 0.0
    9'hx59,9'hx69,        // 0**infinity is 0.0
    9'hx92,               // infinity**neg_infinity is 0.0
    9'hx93,9'hx94: begin  // infinity**-real is 0.0
      power_of_fp = fp_gen_number(4'd6); // return 0
    end
    9'hx35,9'hx36,
    9'hx45,9'hx46,        // -real**0 is 1.0
    9'hx75,9'hx76,
    9'hx85,9'hx86,        // +real**0 is 1.0
    9'hx95,9'hx96: begin  // infinity**0 is 1.0
      power_of_fp = fp_gen_number(4'd8);
    end
    9'hx72,9'hx82,        // +real**neg_infinity
    9'hx79,9'hx89: begin  // +real**infinity
      if(real_l == 1.0)
        power_of_fp = fp_gen_number(4'd8); // return 1.0
      else if(real_l > 1.0) begin
        if(r_fptype == 4'd2)
          power_of_fp = fp_gen_number(4'd6); // return 0.0
        else
          power_of_fp = fp_gen_number(4'd9); // infinity
      end
      else begin
        if(r_fptype == 4'd9)
          power_of_fp = fp_gen_number(4'd6); // return 0.0
        else
          power_of_fp = fp_gen_number(4'd9); // infinity
      end
    end
    9'h133,9'h143,        // -real**-integer
    9'h138,9'h148,        // -real**+integer
    9'hx73,9'hx74,
    9'hx83,9'hx84,        // +real**-real
    9'hx77,9'hx78,
    9'hx87,9'hx88: begin  // +real**+real
      log2_abs_l = $bitstoreal(logarithm(32'd2,{1'b0,l[62:0]}));
      log2_result = real_r * log2_abs_l;
      signed_reg_log2_result = fp64_to_integer_func($realtobits(log2_result));
      integer_log2_result = signed_reg_log2_result;
      if(    (integer_log2_result != log2_result)
          && (integer_log2_result < 0))
        integer_log2_result = integer_log2_result - 1;
      if(    (integer_log2_result <= log2_result)
          && ((integer_log2_result + 1) > log2_result)) begin     
        log2_result_fraction = log2_result - integer_log2_result;
        if(log2_result_fraction != 0.0) begin
          pow_last = 0;
          term = 1.0;
          pow = term;
          term1 = log2_result_fraction * $bitstoreal(natural_log($realtobits(2.0)));
          i=1;
          while(pow_last != pow) begin
            pow_last = pow;
            term = term * term1 / i;
            pow = pow_last + term;
            i = i+1;
          end
        end
        else
          pow = 1.0;
        pow_bits = $realtobits(pow);
        pow_fptype = fp_class(pow_bits);
        if(    (pow_fptype == 4'd8)
            && ((integer_log2_result < 0) || (integer_log2_result <  11'h400))
            && ((integer_log2_result > -60) || ((-integer_log2_result-52) < 11'h3ff))) begin
          if((integer_log2_result < 0) && (-integer_log2_result >= 11'h3ff)) begin
            // make a denormalized result
            fraction = {1'b1,pow_bits[52-1:0],3'b0} >> (-integer_log2_result - 11'h3fe);
            if(&fraction[52-1:-1])
              power_of_fp = {(l[63] & r_is_an_integer & signed_reg_r),11'h001,52'b0};
            else begin
              if(fraction[-1]) begin // rounding
                if(~|fraction[-2:-3]) begin
                  if(fraction[0])
                    fraction = fraction + 4'b1000;
                end
                else   
                  fraction = fraction + 4'b1000;
              end
              power_of_fp = {(l[63] & r_is_an_integer & signed_reg_r),
                             11'h000,fraction[52-1:0]};
            end
          end
          else begin
            signed_exp = integer_log2_result - 1;
            power_of_fp = {(l[63] & r_is_an_integer & signed_reg_r),
                           ~signed_exp[11-1],
                           signed_exp[11-2:0],
                           pow_bits[52-1:0]};
          end
        end
        else if(    (pow_fptype == 4'd8)
                 && ((integer_log2_result > 0) && (integer_log2_result >=  11'h400))) begin
          if(l[63])
            power_of_fp = fp_gen_number(4'd1); // return NAN
          else
            power_of_fp = fp_gen_number(4'd9); // infinity
        end 
        else if(    (pow_fptype == 4'd8)
                 && ((integer_log2_result < 0) && ((-integer_log2_result-52) >= 11'h3ff))) begin
          if(l[63])
            power_of_fp = fp_gen_number(4'd1); // return NAN
          else
            power_of_fp = fp_gen_number(4'd6); // return 0
        end
        else
          power_of_fp = fp_gen_number(4'd1); // return NAN 
      end
      else begin
        if(l[63])
          power_of_fp = fp_gen_number(4'd1); // return NAN
        else if(log2_result > 0)
          power_of_fp = fp_gen_number(4'd9); // infinity
        else
          power_of_fp = fp_gen_number(4'd6); // return 0
      end 
    end
    default: begin
	power_of_fp = fp_gen_number(4'd1); // return NAN
    end
    endcase
  end
  else
    power_of_fp = 64'bx;
end
endfunction // power_of_fp


//  -- NOT Synthisable
function [63:0] power_of;
   input [63:0] l;
   input signed [31:0] r; // integer
   integer integer_r;
   reg [63:0] fp_r_bits;
begin
  if(&(r | ~r) === 1'b1) begin
    integer_r = r;
    fp_r_bits = $realtobits($itor(integer_r));
    power_of = power_of_fp(l,fp_r_bits);
  end
  else
    power_of = 64'bx;
end
endfunction // power_of


//  -- Synthisable
function [63:0] e64;
input	nothing;
begin
    e64 = 64'h4005_bf0a_8b14_5769;
end
endfunction // e64


//  -- Synthisable
function [63:0] pi64;
input	nothing;
begin
    pi64 = 64'h400921fb54442d18;
end
endfunction // pi64


//  -- Angle functions
//  -- sine, cosine, tangent, arc_sine, arc_cosine, arc_tangent
 

//  -- NOT Synthisable
function [63:0] sine;
   input [63:0] angle; // in radians
   real real_angle,real_angle_squared,term,sin,sin_last;
   integer i;
begin
  if(&(angle | ~angle) === 1'b1) begin
    real_angle = $bitstoreal(angle);
    real_angle_squared = real_angle * real_angle;
    sin_last = 0;
    term = real_angle;
    sin = term;
    i=1;
    while(sin_last != sin) begin
        sin_last = sin;
        term = term * real_angle_squared / ((i+1) * (i+2));
        if(i%4 == 1)
            sin = sin_last - term;
	else
	    sin = sin_last + term;
        i = i+2;
    end
    sine = $realtobits(sin);
  end
  else
    sine = 64'bx;
end
endfunction // sine

 
//  -- NOT Synthisable
function [63:0] cosine;
   input [63:0] angle; // in radians
   real real_angle,real_angle_squared,term,cos,cos_last;
   integer i;
begin
  if(&(angle | ~angle) === 1'b1) begin
    real_angle = $bitstoreal(angle);
    real_angle_squared = real_angle * real_angle;
    cos_last = 0;
    term = 1;
    cos = term;
    i=0;
    while(cos_last != cos) begin
        cos_last = cos;
        term = term * real_angle_squared / ((i+1) * (i+2));
        if(i%4 == 0)
            cos = cos_last - term;
	else
	    cos = cos_last + term;
        i = i+2;
    end
    cosine = $realtobits(cos);
  end
  else
    cosine = 64'bx;
end
endfunction // cosine

 
//  -- NOT Synthisable
function [63:0] tangent;
   input [63:0] angle; // in radians
   real tan;
begin
  if(&(angle | ~angle) === 1'b1) begin
    tan = $bitstoreal(sine(angle)) / $bitstoreal(cosine(angle));
    tangent = $realtobits(tan);
  end
  else
    tangent = 64'bx;
end
endfunction // tangent

 
//  -- NOT Synthisable
function [63:0] arc_sine;
   input [63:0] arg; // in radians
   real real_arg,real_arg_squared,term,arc_sin,arc_sin_last;
   integer i;
   reg [63:0] arc_sin_bits;
begin
  if(&(arg | ~arg) === 1'b1) begin
    real_arg = $bitstoreal(arg);
    real_arg_squared = real_arg * real_arg;
    if(real_arg_squared < 1.0) begin
        arc_sin_last = 0;
        term = real_arg;
        arc_sin = term;
        i=1;
        while(arc_sin_last != arc_sin) begin
            arc_sin_last = arc_sin;
            term = i * i * term * real_arg_squared / ((i+1) * (i+2));
	    arc_sin = arc_sin_last + term;
            i = i+2;
        end
        arc_sine = $realtobits(arc_sin);
    end
    else if(real_arg_squared == 1.0) begin
	term = $bitstoreal(pi64(1)) / 2;
	arc_sin_bits = $realtobits(term);
        arc_sine = {arg[63],arc_sin_bits[62:0]};
    end
    else
        arc_sine = fp_gen_number(4'd1);
  end
  else
    arc_sine = 64'bx;
end
endfunction // arc_sine


//  -- NOT Synthisable
function [63:0] arc_cosine;
   input [63:0] arg; // in radians
   real real_arg,real_arg_squared,arc_cos;
begin
  if(&(arg | ~arg) === 1'b1) begin
    real_arg = $bitstoreal(arg);
    real_arg_squared = real_arg * real_arg;
    if(real_arg_squared <= 1.0) begin	
	arc_cos = ($bitstoreal(pi64(1)) / 2) - $bitstoreal(arc_sine(arg));
        arc_cosine = $realtobits(arc_cos);
    end
    else
        arc_cosine = fp_gen_number(4'd1);
  end
  else
    arc_cosine = 64'bx;
end
endfunction // arc_cosine

 
//  -- NOT Synthisable
function [63:0] arc_tangent;
   input [63:0] arg; // in radians
   real real_arg,real_arg_squared,term,term1,arc_tan,arc_tan_last;
   integer i;
   reg [63:0] arc_tan_bits;
begin
  if(&(arg | ~arg) === 1'b1) begin
    real_arg = $bitstoreal(arg);
    real_arg_squared = real_arg * real_arg;
    if($bitstoreal(absolute(arg)) <= 1.0 ) begin
        arc_tan_last = 0;
        term = real_arg / ( 1.0 + real_arg_squared);
        term1 = real_arg_squared / (1.0 + real_arg_squared);
        arc_tan = term;
        i=1;
        while(arc_tan_last != arc_tan) begin
            arc_tan_last = arc_tan;
            term = (i+1) * term * term1 / (i+2);
            arc_tan = arc_tan_last + term;
            i = i+2;
        end
        arc_tangent = $realtobits(arc_tan);
    end
    else begin
        arc_tan_last = 0;
        term = 1.0 / real_arg;
        arc_tan = -term;
        i=1;
        while(arc_tan_last != arc_tan) begin
            arc_tan_last = arc_tan;
            term = i * term / ((i+2) * real_arg_squared);
            if(i%4==1)
                arc_tan = arc_tan_last + term;
            else
                arc_tan = arc_tan_last - term;
            i = i+2;
        end
	if(arg[63])
            arc_tan = arc_tan - ($bitstoreal(pi64(1)) / 2);
        else
            arc_tan = arc_tan + ($bitstoreal(pi64(1)) / 2);
        arc_tangent = $realtobits(arc_tan);
    end
  end
  else
    arc_tangent = 64'bx;
end
endfunction // arc_tangent


//  -- Complex function
//  -- NOT Synthisable
function [2*64-1:0] complex_add;
   input [63:0] l_r, l_i, r_r, r_i;
   reg [63:0] res_r,res_i;
begin
    res_r = add(l_r,r_r);
    res_i = add(l_i,r_i);
    complex_add = {res_r,res_i};
end
endfunction // complex_add

//  -- Complex function
//  -- NOT Synthisable
function [2*64-1:0] complex_subtract;
   input [63:0] l_r, l_i, r_r, r_i;
   reg [63:0] res_r,res_i;
begin
    res_r = subtract(l_r,r_r);
    res_i = subtract(l_i,r_i);
    complex_subtract = {res_r,res_i};
end
endfunction // complex_subtract

//  -- Complex function
//  -- NOT Synthisable
function [2*64-1:0] complex_multiply;
   input [63:0] l_r, l_i, r_r, r_i;
   real real_l_r,real_r_r,real_res_r;
   real real_l_i,real_r_i,real_res_i;
   reg	[63:0] res_r,res_i;
begin
   real_l_r = $bitstoreal(l_r);
   real_l_i = $bitstoreal(l_i);
   real_r_r = $bitstoreal(r_r);
   real_r_i = $bitstoreal(r_i);
   real_res_r = (real_l_r * real_r_r) + (real_l_i * real_r_i);
   real_res_i = (real_l_r * real_r_i) + (real_l_i * real_r_r);
   res_r = $realtobits(real_res_r);
   res_i = $realtobits(real_res_i);
   complex_multiply = {res_r,res_i};
end
endfunction // complex_multiply

//  -- Complex function
//  -- NOT Synthisable
function [2*64-1:0] complex_divide;
   input [63:0] l_r, l_i, r_r, r_i;
   real real_l_r,real_r_r,real_res_r;
   real real_l_i,real_r_i,real_res_i;
   reg	[63:0] res_r,res_i;
begin
   real_l_r = $bitstoreal(l_r);
   real_l_i = $bitstoreal(l_i);
   real_r_r = $bitstoreal(r_r);
   real_r_i = $bitstoreal(r_i);
   real_res_r = ((real_l_r * real_r_r) + (real_l_i * real_r_i))
		/ ((real_l_i * real_l_i) + (real_r_i * real_r_i));
   real_res_i = ((real_l_i * real_r_r) + (real_l_r * real_r_i))
		/ ((real_l_i * real_l_i) + (real_r_i * real_r_i));
   res_r = $realtobits(real_res_r);
   res_i = $realtobits(real_res_i);
   complex_divide = {res_r,res_i};
end
endfunction // complex_divide


//  -- Hyperbolic functions?


//  -- compare functions
//  -- =, !=, >=, <=, <, >
//  -- NOT Synthisable
function equal;
   input [63:0] l, r;
   real real_l,real_r;
begin
   real_l = $bitstoreal(l);
   real_r = $bitstoreal(r);
   equal = (real_l == real_r);
end
endfunction // equal
  
//  -- NOT Synthisable
function not_equal;
   input [63:0] l, r;
   real real_l,real_r;
begin
   real_l = $bitstoreal(l);
   real_r = $bitstoreal(r);
   not_equal = (real_l != real_r);
end
endfunction // not_equal
  
//  -- NOT Synthisable
function greater_than_or_equal;
   input [63:0] l, r;
   real real_l,real_r;
begin
   real_l = $bitstoreal(l);
   real_r = $bitstoreal(r);
   greater_than_or_equal = (real_l >= real_r);
end
endfunction // greater_than_or_equal
  
//  -- NOT Synthisable
function less_than_or_equal;
   input [63:0] l, r;
   real real_l,real_r;
begin
   real_l = $bitstoreal(l);
   real_r = $bitstoreal(r);
   less_than_or_equal = (real_l <= real_r);
end
endfunction // less_than_or_equal
  
//  -- NOT Synthisable
function greater_than;
   input [63:0] l, r;
   real real_l,real_r;
begin
   real_l = $bitstoreal(l);
   real_r = $bitstoreal(r);
   greater_than = (real_l > real_r);
end
endfunction // greater_than
  
//  -- NOT Synthisable
function less_than;
   input [63:0] l, r;
   real real_l,real_r;
begin
   real_l = $bitstoreal(l);
   real_r = $bitstoreal(r);
   less_than = (real_l < real_r);
end
endfunction // less_than

  
//  -- conversion functions

//  -- NOT Synthisable
function [63:0] integer_to_fp64_func;
input signed [31:0] arg;
   real real_res;
begin
   real_res = arg;
   integer_to_fp64_func = $realtobits(real_res);
end
endfunction // integer_to_fp64_func

//  -- NOT Synthisable
function integer fp64_to_integer_func;
input [63:0] arg;
   real real_arg;
begin
   real_arg = $bitstoreal(arg);
   fp64_to_integer_func = $rtoi(real_arg);
end
endfunction // fp64_to_integer_func