Merge FPGA repository back into UHD repository

The FPGA codebase was removed from the UHD repository in 2014 to reduce
the size of the repository. However, over the last half-decade, the
split between the repositories has proven more burdensome than it has
been helpful. By merging the FPGA code back, it will be possible to
create atomic commits that touch both FPGA and UHD codebases. Continuous
integration testing is also simplified by merging the repositories,
because it was previously difficult to automatically derive the correct
UHD branch when testing a feature branch on the FPGA repository.

This commit also updates the license files and paths therein.

We are therefore merging the repositories again. Future development for
FPGA code will happen in the same repository as the UHD host code and
MPM code.

== Original Codebase and Rebasing ==

The original FPGA repository will be hosted for the foreseeable future
at its original local location: https://github.com/EttusResearch/fpga/

It can be used for bisecting, reference, and a more detailed history.

The final commit from said repository to be merged here is
05003794e2da61cabf64dd278c45685a7abad7ec. This commit is tagged as
v4.0.0.0-pre-uhd-merge.

If you have changes in the FPGA repository that you want to rebase onto
the UHD repository, simply run the following commands:

- Create a directory to store patches (this should be an empty
  directory):

    mkdir ~/patches

- Now make sure that your FPGA codebase is based on the same state as
  the code that was merged:

    cd src/fpga # Or wherever your FPGA code is stored
    git rebase v4.0.0.0-pre-uhd-merge

  Note: The rebase command may look slightly different depending on what
  exactly you're trying to rebase.

- Create a patch set for your changes versus v4.0.0.0-pre-uhd-merge:

    git format-patch v4.0.0.0-pre-uhd-merge -o ~/patches

  Note: Make sure that only patches are stored in your output directory.
  It should otherwise be empty. Make sure that you picked the correct
  range of commits, and only commits you wanted to rebase were exported
  as patch files.

- Go to the UHD repository and apply the patches:

    cd src/uhd # Or wherever your UHD repository is stored
    git am --directory fpga ~/patches/*
    rm -rf ~/patches # This is for cleanup

== Contributors ==

The following people have contributed mainly to these files (this list
is not complete):

Co-authored-by: Alex Williams <alex.williams@ni.com>
Co-authored-by: Andrej Rode <andrej.rode@ettus.com>
Co-authored-by: Ashish Chaudhari <ashish@ettus.com>
Co-authored-by: Ben Hilburn <ben.hilburn@ettus.com>
Co-authored-by: Ciro Nishiguchi <ciro.nishiguchi@ni.com>
Co-authored-by: Daniel Jepson <daniel.jepson@ni.com>
Co-authored-by: Derek Kozel <derek.kozel@ettus.com>
Co-authored-by: EJ Kreinar <ej@he360.com>
Co-authored-by: Humberto Jimenez <humberto.jimenez@ni.com>
Co-authored-by: Ian Buckley <ian.buckley@gmail.com>
Co-authored-by: Jörg Hofrichter <joerg.hofrichter@ni.com>
Co-authored-by: Jon Kiser <jon.kiser@ni.com>
Co-authored-by: Josh Blum <josh@joshknows.com>
Co-authored-by: Jonathon Pendlum <jonathan.pendlum@ettus.com>
Co-authored-by: Martin Braun <martin.braun@ettus.com>
Co-authored-by: Matt Ettus <matt@ettus.com>
Co-authored-by: Michael West <michael.west@ettus.com>
Co-authored-by: Moritz Fischer <moritz.fischer@ettus.com>
Co-authored-by: Nick Foster <nick@ettus.com>
Co-authored-by: Nicolas Cuervo <nicolas.cuervo@ettus.com>
Co-authored-by: Paul Butler <paul.butler@ni.com>
Co-authored-by: Paul David <paul.david@ettus.com>
Co-authored-by: Ryan Marlow <ryan.marlow@ettus.com>
Co-authored-by: Sugandha Gupta <sugandha.gupta@ettus.com>
Co-authored-by: Sylvain Munaut <tnt@246tNt.com>
Co-authored-by: Trung Tran <trung.tran@ettus.com>
Co-authored-by: Vidush Vishwanath <vidush.vishwanath@ettus.com>
Co-authored-by: Wade Fife <wade.fife@ettus.com>

5 files changed, 646 insertions, 0 deletions

diff --git a/fpga/usrp3/sim/general/Makefile.srcs b/fpga/usrp3/sim/general/Makefile.srcs
new file mode 100644
index 000000000..90ce1a028
--- /dev/null
+++ b/fpga/usrp3/sim/general/Makefile.srcs
@@ -0,0 +1,13 @@
+#
+# Copyright 2015 Ettus Research LLC
+#
+
+##################################################
+# General Simulation Libraries/Headers
+##################################################
+SIM_GENERAL_SRCS = $(abspath $(addprefix $(BASE_DIR)/../sim/general/, \
+sim_clks_rsts.vh \
+sim_exec_report.vh \
+sim_math.vh \
+sim_file_io.svh \
+))
diff --git a/fpga/usrp3/sim/general/sim_clks_rsts.vh b/fpga/usrp3/sim/general/sim_clks_rsts.vh
new file mode 100644
index 000000000..62de98faa
--- /dev/null
+++ b/fpga/usrp3/sim/general/sim_clks_rsts.vh
@@ -0,0 +1,94 @@
+//
+// Copyright 2015 Ettus Research LLC
+//
+
+// Generates a persistent clock that starts at t=0 and runs forever
+//
+// Usage: `DEFINE_CLK(clk_name,period,duty_cycle)
+// where
+//  - clk_name:   The clock net to be generated
+//  - period:     Period of the clock in simulator ticks
+//  - duty_cycle: Percentage duty cycle
+//
+`define DEFINE_CLK(clk_name, period, duty_cycle) \
+    reg clk_name = 0; \
+    always begin \
+      clk_name = 1; \
+      #((1.0*period)*(1.0*duty_cycle/100)); \
+      clk_name = 0; \
+      #((1.000*period)*(1.0-(1.0*duty_cycle/100))); \
+    end
+
+// Generates a persistent clock that starts at t=0 and runs forever
+//
+// Usage: `DEFINE_CLK(clk_name,period,duty_cycle)
+// where
+//  - clk_name:   The clock net to be generated
+//  - period:     Period of the clock in simulator ticks
+//  - duty_cycle: Percentage duty cycle
+//
+`define DEFINE_DIFF_CLK(clk_name_p, clk_name_n, period, duty_cycle) \
+    reg clk_name_p = 0; \
+    reg clk_name_n = 1; \
+    always begin \
+      clk_name_p = 1; \
+      clk_name_n = 0; \
+      #((1.0*period)*(1.0*duty_cycle/100)); \
+      clk_name_p = 0; \
+      clk_name_n = 1; \
+      #((1.000*period)*(1.0-(1.0*duty_cycle/100))); \
+    end
+
+// Generates a clock that starts at the specified time and runs forever
+//
+// Usage: `DEFINE_LATE_START_CLK(clk_name,period,duty_cycle,start_time,start_time_res)
+// where
+//  - clk_name:   The clock net to be generated
+//  - period:     Period of the clock in simulator ticks
+//  - duty_cycle: Percentage duty cycle
+//  - start_time: Start time for clock in simulator ticks
+//  - start_time_res: Start time resolution (must be > timescale increment and < start_time)
+//
+`define DEFINE_LATE_START_CLK(clk_name, period, duty_cycle, start_time, start_time_res) \
+    reg clk_name = 0; \
+    reg clk_name``_locked = 0; \
+    always begin \
+      while (!clk_name``_locked) #start_time_res; \
+      clk_name = 1; \
+      #((1.0*period)*(1.0*duty_cycle/100)); \
+      clk_name = 0; \
+      #((1.000*period)*(1.0-(1.0*duty_cycle/100))); \
+    end \
+    initial begin \
+        #(start_time) clk_name``_locked = 1; \
+    end
+
+// Generates an active high reset
+//
+// Usage: `DEFINE_RESET(reset_name,reset_time,reset_duration)
+// where
+//  - reset_name:     The reset net to be generated
+//  - reset_time:     Time at which reset will be asserted (i.e. rst=1)
+//  - reset_duration: Duration of reset assertion
+//
+`define DEFINE_RESET(reset_name, reset_time, reset_duration) \
+    reg reset_name = (reset_time==0); \
+    initial begin \
+        #(reset_time) reset_name = 1; \
+        #(reset_time+reset_duration) reset_name = 0; \
+    end
+
+// Generates an active low reset
+//
+// Usage: `DEFINE_RESET_N(reset_name,reset_time,reset_duration)
+// where
+//  - reset_name:     The reset net to be generated
+//  - reset_time:     Time at which reset will be asserted (i.e. rst=0)
+//  - reset_duration: Duration of reset assertion
+//
+`define DEFINE_RESET_N(reset_name, reset_time, reset_duration) \
+    reg reset_name = (reset_time!=0); \
+    initial begin \
+        #(reset_time) reset_name = 0; \
+        #(reset_time+reset_duration) reset_name = 1; \
+    end
diff --git a/fpga/usrp3/sim/general/sim_exec_report.vh b/fpga/usrp3/sim/general/sim_exec_report.vh
new file mode 100644
index 000000000..a6ede6a25
--- /dev/null
+++ b/fpga/usrp3/sim/general/sim_exec_report.vh
@@ -0,0 +1,138 @@
+//
+// Copyright 2015 Ettus Research LLC
+//
+
+// Initializes state for a test bench.
+// This macro *must be* called within the testbench module but 
+// outside the primary initial block
+// Its sets up boilerplate code for:
+// - Logging to console
+// - Test execution tracking
+// - Gathering test results
+// - Bounding execution time based on the SIM_TIMEOUT_US vdef
+
+`ifndef SIM_TIMEOUT_US
+`define SIM_TIMEOUT_US 100000 // Default: 100 ms
+`endif
+
+// Usage: `TEST_BENCH_INIT(test_name,min_tc_run_count,ns_per_tick)
+// where
+//  - tb_name:          Name of the testbench. (Only used during reporting)
+//  - min_tc_run_count: Number of test cases in testbench. (Used to detect stalls and inf-loops)
+//  - ns_per_tick:      The time_unit_base from the timescale declaration
+//
+`define TEST_BENCH_INIT(tb_name, min_tc_run_count, ns_per_tick) \
+  localparam sim_time_increment = 100; \
+  reg tc_running = 0; \
+  reg tc_failed = 0; \
+  reg tc_all_done = 0; \
+  real sim_time = 0.0; \
+  integer tc_run_count = 0; \
+  integer tc_pass_count = 0; \
+  \
+  initial begin : tb_timekeeper \
+    #0; \
+    $timeformat(-9, 0, " ns", 10); \
+    $display("========================================================"); \
+    $display("TESTBENCH STARTED: %s", tb_name); \
+    $display("========================================================"); \
+    if (1000.0*`SIM_TIMEOUT_US < sim_time_increment) begin \
+      $error("Total simulation time less than simulation step size!"); \
+    end \
+    tc_running = 0; \
+    tc_failed = 0; \
+    tc_run_count = 0; \
+    tc_pass_count = 0; \
+    while (~tc_all_done & sim_time < 1000.0*`SIM_TIMEOUT_US) begin \
+       #(sim_time_increment); \
+       sim_time += sim_time_increment; \
+    end \
+    $display("========================================================"); \
+    $display("TESTBENCH FINISHED: %s", tb_name); \
+    $display(" - Time elapsed:   %0t%s", $realtime(), (sim_time >= 1000.0*`SIM_TIMEOUT_US) ? " (Timed out!)" : ""); \
+    $display(" - Tests Expected: %0d", min_tc_run_count); \
+    $display(" - Tests Run:      %0d", tc_run_count); \
+    $display(" - Tests Passed:   %0d", tc_pass_count); \
+    $display("Result: %s", ((tc_run_count>=min_tc_run_count)&&(tc_run_count==tc_pass_count)?"PASSED   ":"FAILED!!!")); \
+    $display("========================================================"); \
+    $finish; \
+  end
+
+// Ends test bench. Place after final test.
+//
+// Usage: `TEST_BENCH_DONE
+`define TEST_BENCH_DONE tc_all_done = 1;
+
+// Indicates the start of a test case
+// This macro *must be* called inside the primary initial block
+//
+// Usage: `TEST_CASE_START(test_name)
+// where
+//  - test_name:        The name of the test.
+//
+`define TEST_CASE_START(test_name) \
+  #0; \
+  tc_running = 1; \
+  tc_failed = 0; \
+  tc_run_count = tc_run_count + 1; \
+  $display("[TEST CASE %3d] (t=%09d) BEGIN: %s...", tc_run_count, $time, test_name);
+
+// Indicates the end of a test case
+// This macro *must be* called inside the primary initial block
+// The pass/fail status of test case is determined based on the
+// the user specified outcome and the number of fatal or error
+// ASSERTs triggered in the test case.
+//
+// Usage: `TEST_CASE_DONE(test_result)
+// where
+//  - test_result:  User specified outcome
+//
+`define TEST_CASE_DONE(result) \
+  #0; \
+  tc_running = 0; \
+  $display("[TEST CASE %3d] (t=%09d) DONE... %s", tc_run_count, $time, ((((result)===1'b1)&~tc_failed)?"Passed":"FAILED")); \
+  if (((result)===1'b1)&~tc_failed) tc_pass_count = tc_pass_count + 1;
+
+// Wrapper around a an assert.
+// ASSERT_FATAL throws an error assertion and halts the simulator
+// if cond is not satisfied
+//
+// Usage: `ASSERT_FATAL(cond,msg)
+// where
+//  - cond: Condition for the assert
+//  - msg:  Message for the assert
+//
+`define ASSERT_FATAL(cond, msg) \
+  assert(cond) else begin \
+    tc_failed = 1; \
+    $error(msg); \
+    tc_all_done = 1; \
+    #(sim_time_increment); \
+  end
+
+// Wrapper around a an assert.
+// ASSERT_ERROR throws an error assertion and fails the test case
+// if cond is not satisfied. The simulator will *not* halt
+//
+// Usage: `ASSERT_ERROR(cond,msg)
+// where
+//  - cond: Condition for the assert
+//  - msg:  Message for the assert
+//
+`define ASSERT_ERROR(cond, msg) \
+  assert(cond) else begin \
+    tc_failed = 1; \
+    $error(msg); \
+  end
+
+// Wrapper around a an assert.
+// ASSERT_WARNING throws an warning assertion but does not fail the
+// test case if cond is not satisfied. The simulator will *not* halt
+//
+// Usage: `ASSERT_WARNING(cond,msg)
+// where
+//  - cond: Condition for the assert
+//  - msg:  Message for the assert
+//
+`define ASSERT_WARN(cond, msg) \
+  assert(cond) else $warning(msg);
diff --git a/fpga/usrp3/sim/general/sim_file_io.svh b/fpga/usrp3/sim/general/sim_file_io.svh
new file mode 100644
index 000000000..97fadd741
--- /dev/null
+++ b/fpga/usrp3/sim/general/sim_file_io.svh
@@ -0,0 +1,125 @@
+//
+// Copyright 2015 Ettus Research LLC
+//
+
+`ifndef WORKING_DIR
+  `define WORKING_DIR "."
+`endif
+
+`define ABSPATH(name) {`WORKING_DIR, "/", name}
+
+typedef enum {
+  READ, WRITE, APPEND
+} fopen_mode_t;
+
+typedef enum {
+  HEX, DEC, OCT, BIN, FLOAT
+} fformat_t;
+
+// Create a handle to a data_file with
+// - FILENAME: Name of the file
+// - FORMAT: Data format (HEX, DEC, OCT, BIN, FLOAT)
+// - DWIDTH: Width of each element stored in the file (one line per word)
+//
+
+//TODO: We would ideally use a class but that is not
+//      supported by most simulators.
+interface data_file_t #(
+  parameter FILENAME = "test.hex", 
+  parameter FORMAT = HEX, 
+  parameter DWIDTH = 64
+) (input clk);
+  bit     is_open;
+  integer handle;
+
+  // Open the data file for reading or writing.
+  //
+  // Usage: open(mode)
+  // where
+  //  - mode: RW mode (Choose from: READ, WRITE, APPEND)
+  //
+  function open(fopen_mode_t mode = READ);
+    if (mode == APPEND)
+      handle  = $fopen(`ABSPATH(FILENAME), "a");
+    else if (mode == WRITE)
+      handle  = $fopen(`ABSPATH(FILENAME), "w");
+    else
+      handle  = $fopen(`ABSPATH(FILENAME), "r");
+
+    if (handle == 0) begin
+      $error("Could not open file: %s", `ABSPATH(FILENAME));
+      $finish();
+    end
+    is_open = 1;
+  endfunction
+
+
+  // Close an open data file. No-op if file isn't already open
+  //
+  // Usage: close()
+  //
+  function close();
+    $fclose(handle);
+    handle  = 0;
+    is_open = 0;
+  endfunction
+
+  // Is end-of-file reached.
+  //
+  // Usage: is_eof() Returns eof
+  // where
+  // - eof: A boolean
+  //
+  function logic is_eof();
+    return ($feof(handle));
+  endfunction
+
+  // Read a line from the datafile
+  //
+  // Usage: readline() Returns data
+  // where
+  // - data: A logic array of width DWIDTH containing the read word
+  //
+  function logic [DWIDTH-1:0] readline();
+    automatic logic [DWIDTH-1:0] word = 64'h0;
+    automatic integer status;
+
+    if (FORMAT == HEX)
+      status = $fscanf(handle, "%x\n", word);
+    else if (FORMAT == DEC)
+      status = $fscanf(handle, "%d\n", word);
+    else if (FORMAT == OCT)
+      status = $fscanf(handle, "%o\n", word);
+    else if (FORMAT == BIN)
+      status = $fscanf(handle, "%b\n", word);
+    else if (FORMAT == DEC)
+      status = $fscanf(handle, "%g\n", word);
+    else
+      $error("Invalid format");
+
+    return word;
+  endfunction
+
+  // Write a line to the datafile
+  //
+  // Usage: writeline(data) 
+  // where
+  // - data: A logic array of width DWIDTH to write to the file
+  //
+  function void writeline(logic [DWIDTH-1:0] word);
+    if (FORMAT == HEX)
+      $fdisplay(handle, "%x", word);
+    else if (FORMAT == DEC)
+      $fdisplay(handle, "%d", word);
+    else if (FORMAT == OCT)
+      $fdisplay(handle, "%o", word);
+    else if (FORMAT == BIN)
+      $fdisplay(handle, "%b", word);
+    else if (FORMAT == DEC)
+      $fdisplay(handle, "%g", word);
+    else
+      $error("Invalid format");
+  endfunction
+
+endinterface
+
diff --git a/fpga/usrp3/sim/general/sim_math.vh b/fpga/usrp3/sim/general/sim_math.vh
new file mode 100644
index 000000000..a08870740
--- /dev/null
+++ b/fpga/usrp3/sim/general/sim_math.vh
@@ -0,0 +1,276 @@
+// All code take from the HDLCon paper:
+// "Verilog Transcendental Functions for Numerical Testbenches"
+//
+// Authored by:
+// Mark G. Arnold marnold@co.umist.ac.uk,
+// Colin Walter c.walter@co.umist.ac.uk
+// Freddy Engineer freddy.engineer@xilinx.com
+//
+
+
+
+// The sine function is approximated with a polynomial which works
+// for -π/2 < x < π/2. (This polynomial, by itself, was used as a
+// Verilog example in [2]; unfortunately there was a typo with the
+// coefficients. The correct coefficients together with an error
+// analysis are given in [3].)   For arguments outside of -π/2 < x < π/2,
+// the identities sin(x) = -sin(-x) and sin(x) = -sin(x-π) allow the
+// argument to be shifted to be within this range.  The latter identity
+// can be applied repeatedly.  Doing so could cause inaccuracies for
+// very large arguments, but in practice the errors are acceptable
+// if the Verilog simulator uses double-precision floating point.
+
+function real sin;
+  input x;
+  real x;
+  real x1,y,y2,y3,y5,y7,sum,sign;
+  begin
+    sign = 1.0;
+    x1 = x;
+    if (x1<0)
+      begin
+        x1 = -x1;
+        sign = -1.0;
+      end
+    while (x1 > 3.14159265/2.0)
+      begin
+        x1 = x1 - 3.14159265;
+        sign = -1.0*sign;
+      end
+    y = x1*2/3.14159265;
+    y2 = y*y;
+    y3 = y*y2;
+    y5 = y3*y2;
+    y7 = y5*y2;
+    sum = 1.570794*y - 0.645962*y3 +
+           0.079692*y5 - 0.004681712*y7;
+    sin = sign*sum;
+  end
+endfunction
+
+// The cosine and tangent are computed from the sine:
+function real cos;
+  input x;
+  real x;
+  begin
+    cos = sin(x + 3.14159265/2.0);
+  end
+endfunction
+
+
+function real tan;
+  input x;
+  real x;
+  begin
+    tan = sin(x)/cos(x);
+  end
+endfunction
+
+// The base-two exponential (antilogarithm) function, 2x, is computed by
+// examining the bits of the argument, and for those bits of the argument
+// that are 1, multiplying the result by the corresponding power of a base
+//  very close to one.  For example,  if there were only two bits after
+// the radix point, the base would be the fourth root of two, 1.1892.
+// This number is squared on each iteration:  1.4142,  2.0,  4.0,  16.0.
+// So, if x is 101.112, the function computes 25.75 as 1.1892*1.4142*2.0*16.0 = 53.81.
+// In general, for k bits of precision, the base would be the 2k root of two.
+// Since we need about 23 bits of accuracy for our function, the base we use
+// is the 223 root of two, 1.000000082629586.  This constant poses a problem
+// to some Verilog parsers, so we construct it in two parts.  The following
+// function computes the appropriate root of two by repeatedly squaring this constant:
+
+function real rootof2;
+  input n;
+  integer n;
+  real power;
+  integer i;
+
+  begin
+    power = 0.82629586;
+    power = power / 10000000.0;
+    power = power + 1.0;
+    i = -23;
+
+    if (n >= 1)
+      begin
+        power = 2.0;
+        i = 0;
+      end
+
+    for (i=i; i< n; i=i+1)
+      begin
+        power = power * power;
+      end
+    rootof2 = power;
+  end
+endfunction // if
+
+// This function is used for computing both antilogarithms and logarithms.
+// This routine is never called with n less than -23, thus no validity check
+// need be performed. When n>0, the exponentiation begins with 2.0 in order to
+// improve accuracy.
+// For computing the antilogarithm, we make use of the identity ex = 2x/ln(2),
+// and then proceed as in the example above.  The constant 1/ln(2) = 1.44269504.
+// Here is the natural exponential function:
+
+function real exp;
+  input x;
+  real x;
+  real x1,power,prod;
+  integer i;
+  begin
+    x1 = fabs(x)*1.44269504;
+    if (x1 > 255.0)
+      begin
+        exp = 0.0;
+        if (x>0.0)
+          begin
+            $display("exp illegal argument:",x);
+            $stop;
+          end
+      end
+    else
+      begin
+        prod = 1.0;
+        power = 128.0;
+        for (i=7; i>=-23; i=i-1)
+          begin
+            if (x1 > power)
+              begin
+                prod = prod * rootof2(i);
+                x1 = x1 - power;
+              end
+            power = power / 2.0;
+          end
+        if (x < 0)
+          exp = 1.0/prod;
+        else
+          exp = prod;
+      end
+  end
+endfunction // fabs
+
+// The function prints an error message if the argument is too large
+// (greater than about 180).  All error messages in this package are
+// followed by $stop  to allow the designer to use the debugging
+// features of Verilog to determine the cause of the error, and
+// possibly to resume the simulation.  An argument of less than
+// about –180 simply returns zero with no error.  The main loop
+// assumes a positive argument.  A negative argument is computed as 1/e-x.
+// The logarithm function prints an error message for arguments less
+// than or equal to zero because the real-valued logarithm is not
+// defined for such arguments.  The loop here requires an argument
+// greater than or equal to one.  For arguments between zero and one,
+// this code uses the identity ln(1/x) = -ln(x).
+
+function real log;
+  input x;
+  real x;
+  real re,log2;
+  integer i;
+  begin
+    if (x <= 0.0)
+      begin
+        $display("log illegal argument:",x);
+        $stop;
+        log = 0;
+      end
+    else
+      begin
+        if (x<1.0)
+          re = 1.0/x;
+        else
+          re = x;
+        log2 = 0.0;
+        for (i=7; i>=-23; i=i-1)
+          begin
+            if (re > rootof2(i))
+              begin
+                re = re/rootof2(i);
+                log2 = 2.0*log2 + 1.0;
+              end
+            else
+              log2 = log2*2;
+          end
+        if (x < 1.0)
+          log = -log2/12102203.16;
+        else
+          log = log2/12102203.16;
+      end
+    end
+endfunction
+
+// The code only divides re by rootof2(i) when the re is larger
+// (so that the quotient will be greater than 1.0). Each time
+// such a division occurs, a bit that is 1 is recorded in the
+// whole number result (multiply by 2 and add 1).  Otherwise,
+// a zero is recorded (multiply by 2).  At the end of the loop,
+// log2 will contain 223 log2|x|.  We divide by 223 and use the
+// identity ln(x) = log2(x)/log2(e).  The constant 12102203.16 is 223  log2(e).
+// The log(x) and exp(x)functions are used to implement the pow(x,y) and sqrt(x) functions:
+
+function real pow;
+  input x,y;
+  real x,y;
+  begin
+    if (x<0.0)
+      begin
+        $display("pow illegal argument:",x);
+        $stop;
+      end
+    pow = exp(y*log(x));
+  end
+endfunction
+
+function real sqrt;
+  input x;
+  real x;
+  begin
+    if (x<0.0)
+      begin
+        $display("sqrt illegal argument:",x);
+        $stop;
+      end
+    sqrt = exp(0.5*log(x));
+  end
+endfunction
+
+// The arctangent [3,7] is computed as a continued fraction,
+// using the identities tan-1(x) = -tan-1(-x) and tan-1(x) = π/2 - tan-1(1/x)
+// to reduce the range to 0 < x < 1:
+
+function real atan;
+  input x;
+  real x;
+  real x1,x2,sign,bias;
+  real d3,s3;
+  begin
+    sign = 1.0;
+    bias = 0.0;
+    x1 = x;
+    if (x1 < 0.0)
+      begin
+        x1 = -x1;
+        sign = -1.0;
+      end
+    if (x1 > 1.0)
+      begin
+        x1 = 1.0/x1;
+        bias = sign*3.14159265/2.0;
+        sign = -1.0*sign;
+      end
+    x2 = x1*x1;
+    d3 = x2 + 1.44863154;
+    d3 = 0.26476862 / d3;
+    s3 = x2 + 3.3163354;
+    d3 = s3 - d3;
+    d3 = 7.10676 / d3;
+    s3 = 6.762139 + x2;
+    d3 = s3 - d3;
+    d3 = 3.7092563 / d3;
+    d3 = d3 + 0.17465544;
+    atan = sign*x1*d3+bias;
+  end
+endfunction
+
+// The other functions (asin(x) and acos(x)) are computed from the arctangent.