#! /usr/bin/python
###
### Generate exhaustive tests for floating-point conversions.
###
### (c) 2024 Straylight/Edgeware
###

###----- Licensing notice ---------------------------------------------------
###
### This file is part of the mLib utilities library.
###
### mLib is free software: you can redistribute it and/or modify it under
### the terms of the GNU Library General Public License as published by
### the Free Software Foundation; either version 2 of the License, or (at
### your option) any later version.
###
### mLib is distributed in the hope that it will be useful, but WITHOUT
### ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
### FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
### License for more details.
###
### You should have received a copy of the GNU Library General Public
### License along with mLib.  If not, write to the Free Software
### Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
### USA.

###--------------------------------------------------------------------------
### Imports.

import sys as SYS
import optparse as OP
import random as R
if SYS.version_info >= (3,):
  from io import StringIO
  xrange = range
  def iterkeys(d): return d.keys()
else:
  from cStringIO import StringIO
  def iterkeys(d): return d.iterkeys()

###--------------------------------------------------------------------------
### Utilities.

def bit(k): "Return an integer with just bit K set."; return 1 << k
def mask(k): "Return an integer with bits 0 to K - 1 set."; return bit(k) - 1
M32 = mask(32)

def explore(wd, lobits, hibits):
  """
  Return an iterator over various WD-bit values.

  Suppose that a test wants to explore various WD-bit fields, but WD might be
  too large to do this exhaustively.  We assume (reasonably, in the case at
  hand of floating-point formats) that the really interesting bits are those
  at the low and high ends of the field, and test small subfields at the ends
  exhaustively, filling in the bits in the middle with zeros, ones, or random
  data.

  So, the generator behaves as follows.  If WD <= LOBITS + HIBITS + 1 then
  the iterator will yield all WD-bit values exhaustively.  Otherwise, it
  yields a sequence which includes all combinations of: every LOBITS-bit
  pattern in the least significant bits; every HIBITS-bit pattern in the most
  significant bits; and all-bits-clear, all-bits-set, and a random pattern in
  the bits in between.
  """

  if wd <= hibits + lobits + 1:
    for i in xrange(bit(wd)): yield i
  else:
    midbit = bit(wd - hibits - lobits)
    hishift = wd - hibits
    m = (midbit - 1) << lobits
    for hi in xrange(bit(hibits)):
      top = hi << hishift
      for lo in xrange(bit(lobits)):
        while True:
          fill = R.randrange(midbit)
          if fill != 0 and fill != midbit - 1: break
        base = lo | top
        yield base
        yield base | (fill << lobits)
        yield base | m

class ExploreParameters (object):
  """
  Simple structure for exploration parameters; see `explore' for background.

  The `explo' and `exphi' attributes are the low and high subfield sizes
  for exponent fields, and `siglo' and `sighi' are the low and high subfield
  sizes for significand fields.
  """
  def __init__(me, explo = 0, exphi = 2, siglo = 1, sighi = 3):
    me.explo, me.exphi = explo, exphi
    me.siglo, me.sighi = siglo, sighi

FMTMAP = {} # maps format names to classes

def with_metaclass(meta, *supers):
  """
  Return an arbitrary instance of the metaclass META.

  The class will have SUPERS (default just `object') as its superclasses.
  This is intended to be used in direct-superclass lists, as a compatibility
  hack, because the Python 2 and 3 syntaxes are wildly different.
  """
  return meta("#<anonymous base %s>" % meta.__name__,
              supers or (object,), dict())

class FormatClass (type):
  """
  Metaclass for format classes.

  If the class defines a `NAME' attribute then register the class in
  `FMTMAP'.
  """
  def __new__(cls, name, supers, dict):
    c = type.__new__(cls, name, supers, dict)
    try: FMTMAP[c.NAME] = c
    except AttributeError: pass
    return c

class IEEEFormat (with_metaclass(FormatClass)):
  """
  Floating point format class.

  Concrete subclasses must define the following class attributes.

    * `HIDDENP' -- true if the format uses a `hidden bit' convention for
      normal numbers.

    * `EXPWD' -- exponent field width, in bits.

    * `PREC' -- precision, in bits, /including/ the hidden bit if any.

  Many useful quantities are derived.

    * `_expbias' is the exponent bias.

    * `_minexp' and `_maxexp' are the minimum and maximum representable
      exponents.

    * `_sigwd' is the width of the significand field.

    * `_paywords' is the number of words required to represent a NaN payload.

    * `_nbits' is the total number of bits in an encoded value.

    * `_rawbytes' is the number of bytes required for an encoded value.
  """

  def __init__(me):
    """
    Initialize an instance.
    """
    me._expbias = mask(me.EXPWD - 1)
    me._maxexp = me._expbias
    me._minexp = 1 - me._expbias
    if me.HIDDENP: me._sigwd = me.PREC - 1
    else: me._sigwd = me.PREC

    me._paywords = (me._sigwd + 29)//32

    me._nbits = 1 + me.EXPWD + me._sigwd
    me._rawbytes = (me._nbits + 7)//8

  def decode(me, x):
    """
    Decode the encoded floating-point value X, represented as an integer.

    Return five quantities (FLAGS, EXP, FW, FRAC, ERR), corresponding mostly
    to the `struct floatbits' representation, characterizing the value
    encoded in X.

      * FLAGS is a list of flag tokens:

         -- `NEG' if the value is negative;
         -- `ZERO' if the value is exactly zero;
         -- `INF' if the value is infinite;
         -- `SNAN' if the value is a signalling NaN; and/or
         -- `QNAN' if the value is a quiet NaN.

        FLAGS will be empty if the value is a strictly positive finite
        number.

      * EXP is the exponent, as a signed integer.  This will be `None' if the
        value is zero, infinite, or a NaN.

      * FW is the length of the fraction, in 32-bit words.  This will be
        `None' if the value is zero or infinite.

      * FRAC is the fraction or payload.  This will be `None' if the value is
        zero or infinite; otherwise it will be an integer, 0 <= FRAC <
        2^{32FW}.  If the value is a NaN, then the FRAC represents the
        payload, /not/ including the quiet bit, left aligned.  Otherwise,
        FRAC is normalized so that 2^{32FW-1} <= FRAC < 2^{32FW}, and the
        value represented is S FRAC 2^{EXP-32FW}, where S = -1 if `NEG' is in
        FLAGS, or +1 otherwise.  The represented value is unchanged by
        multiplying or dividing FRAC by an exact power of 2^{32} and
        (respectively) incrementing or decrementing FW to match, but this
        will affect the output data in a way that affects the tests.

      * ERR is a list of error tokens:

          -- `INVAL' if the encoded value is erroneous (though decoding
             continues anyway).

        ERR will be empty if no error occurred.
    """

    ## Extract fields.
    sig = x&mask(me._sigwd)
    biasedexp = (x >> me._sigwd)&mask(me.EXPWD)
    signbit = (x >> (me._sigwd + me.EXPWD))&1
    if not me.HIDDENP: unitbit = sig >> me.PREC - 1

    ## Initialize flag lists.
    flags = []
    err = []

    ## Capture the sign.  This is always relevant.
    if signbit: flags.append("NEG")

    ## If the exponent field is all-bits-set then we have infinity or NaN.
    if biasedexp == mask(me.EXPWD):

      ## If there's no hidden bit then the unit bit should be /set/, but is
      ## /not/ part of the NaN payload -- or even significant for
      ## distinguishing a NaN from an infinity.  If it's clear, signal an
      ## error; if it's set, then clear it so that we don't have to think
      ## about it again.
      if not me.HIDDENP:
        if unitbit: sig &= mask(me._sigwd - 1)
        else: err.append("INVAL")

      ## If the significand is (now) zero, we have an infinity and there's
      ## nothing else to do.
      if not sig:
        flags.append("INF")
        frac = fw = exp = None

      ## Otherwise determine the NaN flavour and extract the payload.
      else:
        if sig&bit(me.PREC - 2): flags.append("QNAN")
        else: flags.append("SNAN")
        shift = 32*me._paywords + 2 - me.PREC
        frac = (sig&mask(me.PREC - 2)) << shift
        exp = None
        fw = me._paywords

    ## Otherwise we have a finite number.  We handle all of these together.
    else:

      ## If there's no hidden bit, then check that the unit bit matches the
      ## exponent: it should be clear if the exponent field is all-bits-zero
      ## (zero or subnormal numbers), and set otherwise (normal numbers).  If
      ## this isn't the case, signal an error, but continue.  We'll normalize
      ## the number correctly as we go.
      if not me.HIDDENP:
        if (not biasedexp and unitbit) or (biasedexp and not unitbit):
          err.append("INVAL")

      ## If the exponent is all-bits-zero then set it to 1; otherwise, if the
      ## format uses a hidden bit then force the unit bit of our significand
      ## on.  The absolute value is now exactly
      ##
      ##        2^{biasedexp-_expbias-PREC+1} sig
      ##
      ## in all cases.
      if not biasedexp: biasedexp = 1
      elif me.HIDDENP: sig |= bit(me._sigwd)

      ## If the significand is now zero then the value must be zero.
      if not sig:
        flags.append("ZERO")
        frac = fw = exp = None

      ## Otherwise we have a nonzero finite value, which might need
      ## normalization.
      else:
        sigwd = sig.bit_length()
        fw = (sigwd + 31)//32
        exp = biasedexp - me._expbias - me.PREC + sigwd + 1
        frac = sig << (32*fw - sigwd)

    ## All done.
    return flags, exp, frac, fw, err

  def _dump_as_bytes(me, var, x, wd):
    """
    Dump an assignment to VAR of X as a WD-byte binary string.

    Print, on standard output, an assignment `VAR = ...' giving the value of
    X, in hexadecimal, split with spaces into groups of 8 digits from the
    right.
    """

    if not wd:
      print("%s = #empty" % var)
    else:
      out = StringIO()
      for i in xrange(wd - 1, -1, -1):
        out.write("%02x" % ((x >> 8*i)&0xff))
        if i and not i%4: out.write(" ")
      print("%s = %s" % (var, out.getvalue()))

  def _dump_flags(me, var, flags, zero = "0"):
    """
    Dump an assignment to VAR of FLAGS as a list of flags.

    Print, on standard output, an assignment `VAR = ...' giving the named
    flags.  Print ZERO (default `0') if FLAGS is empty.
    """

    if flags: print("%s = %s" % (var, " | ".join(flags)))
    else: print("%s = %s" % (var, zero))

  def genenc(me, ep = ExploreParameters()):
    """
    Print, on standard output, tests of encoding floating-point values.

    The tests will cover positive and negative values, with the exponent and
    signficand fields explored according to the parameters EP.
    """

    print("[enc%s]" % me.NAME)
    for s in xrange(2):
      for e in explore(me.EXPWD, ep.explo, ep.exphi):
        for m in explore(me.PREC - 1, ep.siglo, ep.sighi):
          if not me.HIDDENP and e: m |= bit(me.PREC - 1)
          x = (s << (me.EXPWD + me._sigwd)) | (e << me._sigwd) | m
          flags, exp, frac, fw, err = me.decode(x)
          print("")
          me._dump_flags("f", flags)
          if exp is not None: print("e = %d" % exp)
          if frac is not None:
            while not frac&M32 and fw: frac >>= 32; fw -= 1
            me._dump_as_bytes("m", frac, 4*fw)
          me._dump_as_bytes("z", x, me._rawbytes)
          if err: me._dump_flags("err", err, "OK")

  def gendec(me, ep = ExploreParameters()):
    """
    Print, on standard output, tests of decoding floating-point values.

    The tests will cover positive and negative values, with the exponent and
    signficand fields explored according to the parameters EP.
    """

    print("[dec%s]" % me.NAME)
    for s in xrange(2):
      for e in explore(me.EXPWD, ep.explo, ep.exphi):
        for m in explore(me._sigwd, ep.siglo, ep.sighi):
          x = (s << (me.EXPWD + me._sigwd)) | (e << me._sigwd) | m
          flags, exp, frac, fw, err = me.decode(x)
          print("")
          me._dump_as_bytes("x", x, me._rawbytes)
          me._dump_flags("f", flags)
          if exp is not None: print("e = %d" % exp)
          if frac is not None: me._dump_as_bytes("m", frac, 4*fw)
          if err: me._dump_flags("err", err, "OK")

class MiniFloat (IEEEFormat):
  NAME = "mini"
  EXPWD = 4
  PREC = 4
  HIDDENP = True

class BFloat16 (IEEEFormat):
  NAME = "bf16"
  EXPWD = 8
  PREC = 8
  HIDDENP = True

class Binary16 (IEEEFormat):
  NAME = "f16"
  EXPWD = 5
  PREC = 11
  HIDDENP = True

class Binary32 (IEEEFormat):
  NAME = "f32"
  EXPWD = 8
  PREC = 24
  HIDDENP = True

class Binary64 (IEEEFormat):
  NAME = "f64"
  EXPWD = 11
  PREC = 53
  HIDDENP = True

class Binary128 (IEEEFormat):
  NAME = "f128"
  EXPWD = 15
  PREC = 113
  HIDDENP = True

class DoubleExtended80 (IEEEFormat):
  NAME = "idblext80"
  EXPWD = 15
  PREC = 64
  HIDDENP = False

###--------------------------------------------------------------------------
### Main program.

op = OP.OptionParser \
  (description = "Generate test data for IEEE format encoding and decoding",
   usage = "usage: %prog [-E LO/HI] [-M LO/HI] [[enc|dec]FORMAT]")
for shortopt, longopt, kw in \
  [("-E", "--explore-exponent",
    dict(action = "store", metavar = "LO/HI", dest = "expparam",
         help = "exponent exploration parameters")),
   ("-M", "--explore-significand",
    dict(action = "store", metavar = "LO/HI", dest = "sigparam",
         help = "significand exploration parameters"))]:
  op.add_option(shortopt, longopt, **kw)
opts, args = op.parse_args()

ep = ExploreParameters()
for optattr, loattr, hiattr in [("expparam", "explo", "exphi"),
                                ("sigparam", "siglo", "sighi")]:
  opt = getattr(opts, optattr)
  if opt is not None:
    ok = False
    try: sl = opt.index("/")
    except ValueError: pass
    else:
      try: lo, hi = map(int, (opt[:sl], opt[sl + 1:]))
      except ValueError: pass
      else:
        setattr(ep, loattr, lo)
        setattr(ep, hiattr, hi)
        ok = True
    if not ok: op.error("bad exploration parameter `%s'" % opt)

if not args:
  for fmt in iterkeys(FMTMAP):
    args.append("enc" + fmt)
    args.append("dec" + fmt)
firstp = True
for arg in args:
  tail = fmt = None
  if arg.startswith("enc"): tail = arg[3:]; gen = lambda f: f.genenc(ep)
  elif arg.startswith("dec"): tail = arg[3:]; gen = lambda f: f.gendec(ep)
  if tail is not None: fmt = FMTMAP.get(tail)
  if not fmt: op.error("unknown test group `%s'" % arg)
  if firstp: firstp = False
  else: print("")
  gen(fmt())

###----- That's all, folks --------------------------------------------------