// Computer numbers: Integers and IEEE floating point in C#
// Peter Sestoft * sestoft@itu.dk * 2009-02-26, 2015-11-13, 2019-06-11

// Should add "unit tests" of all System.Math functions, following the
// LaTeX tables.  Challenge: Cannot check for NaNs using ==, so need
// to test as 
// result==expected || Double.IsNaN(result)&&Double.IsNaN(expected)

// Compile with 
//   csc Numbers.cs

using System;
using System.Text;

class Numbers {
  public static void Main(String[] args) {
    SystemInfo();
    Console.WriteLine(7.0%2.0);
    Console.WriteLine(Math.IEEERemainder(7.0,2.0));
    byteExamples();
    intExamples();
    floatExamples();
    doubleExamples();
    ExpFunction();
    // Works only because IEEE handles infinities correctly:
    Console.WriteLine("R(1) = " + R(1));
    Console.WriteLine("R(2) = " + R(2));
    Console.WriteLine("R(3) = " + R(3));
    Console.WriteLine("R(4) = " + R(4));
    // Zeroes
    double zp = 0.0, zn = -0.0;
    Console.WriteLine("0.0 = {0} and -0.0 = {1}", zp, zn);
    Console.WriteLine(zp == zn && 1/zp != 1/zn);  // true
    Console.WriteLine(zp.Equals(zn));
    Console.WriteLine(zp.GetHashCode() + " " + zn.GetHashCode());
    Console.WriteLine(-0.0 < 0.0);
    // Loss of significant digits
    double v = 9876543210.2, w = 9876543210.1;
    Console.WriteLine(v+w-w != v);		  // true
    AllDoubleOperations();
    Point[] ps1 = new Point[] { 
    new Point(2.1, 5.2), new Point(2.2, 5.4), new Point(2.4, 5.8) };
    linear1(ps1);
    linear2(ps1);
    Point[] ps2 = Point.moveAll(ps1, 1E7, 1E7);
    linear1(ps2);
    linear2(ps2);
    Point[] ps3 = Point.moveAll(ps1, 5E7, 5E7);
    linear1(ps3);
    linear2(ps3);
    // Testing sum
    int N = 10000;
    double[] xs = new double[2*N];
    for (int i=0; i<N; i++) {
      xs[2*i] = 1E12;
      xs[2*i+1] = -1;
    }
    double result = N*(-1+1E12);
    Console.WriteLine("\nSum Kahan error = {0}; naive error = {1}",
  		      sumKahan(xs)-result, sum(xs) - result);
    // NaNs; operations preserve the NaN code
    double d = MakeNaN(17);
    displayDouble(d);
    Console.WriteLine(ErrorCode(d));
    Console.WriteLine(ErrorCode(0.0 * d));
    Console.WriteLine(ErrorCode(Math.Log(d)));
    Console.WriteLine(ErrorCode(Math.Sin(d)));
    Console.WriteLine(ErrorCode(Math.Cos(d)));
    Console.WriteLine(ErrorCode(Math.Tan(d)));
    Console.WriteLine(ErrorCode(Math.Exp(d)));
    // Console.WriteLine(ErrorCode(Math.Pow(d,1)));
    // Console.WriteLine(ErrorCode(Math.Pow(1,d)));
    // Console.WriteLine(Math.Pow(d,1));
    TestNanPayload();
  }

  private static void SystemInfo() {
    Console.WriteLine("# OS          {0}", 
      Environment.OSVersion.VersionString);
    Console.WriteLine("# .NET vers.  {0}",   
      Environment.Version);
    Console.WriteLine("# 64-bit OS   {0}",   
      Environment.Is64BitOperatingSystem);
    Console.WriteLine("# 64-bit proc {0}",   
      Environment.Is64BitProcess);
    Console.WriteLine("# CPU         {0}; {1} \"cores\"",   
      Environment.GetEnvironmentVariable("PROCESSOR_IDENTIFIER"),
      Environment.ProcessorCount); 
    Console.WriteLine("# Date        {0:s}", 
      DateTime.Now);
  }

  private static bool NanPayload(Func<double, double> f, String name, int i) { 
    bool ok;
    if (!(ok = ErrorCode(f(MakeNaN(i))) == i))
	Console.WriteLine("Wrong NaN payload in {0}(NaN) at errorbits={1}",
			  name, i);
    return ok;
  }

  private static bool NanPayload(Func<double, double, double> f, 
				 String name, int i) 
  { 
    double d = MakeNaN(i);
    bool ok1, ok2;
    if (!(ok1 = ErrorCode(f(d,2)) == i))
      Console.WriteLine("Wrong NaN payload in {0}(NaN,2) at errorbits={1}",
			name, i);
    if (!(ok2 = ErrorCode(f(2,d)) == i))
      Console.WriteLine("Wrong NaN payload in {0}(2,NaN) at errorbits={1}",
			name, i);
    return ok1 && ok2;
  }

  private static void NanPayload(Func<double, double> f, String name) { 
    for (int i=0; i<1025; i++)
      if (!NanPayload(f, name, i) || !NanPayload(f, name, -i))
	return;
  }

  private static void NanPayload(Func<double, double, double> f, String name) 
  { 
    for (int i=0; i<1025; i++) 
      if (!NanPayload(f, name, i) || !NanPayload(f, name, -i))
	return;
  }

  private static void TestNanPayload() {
    NanPayload(delegate(double x) { return x; }, "identity");
    NanPayload(Math.Abs, "Abs");
    NanPayload(Math.Acos, "Acos");
    NanPayload(Math.Asin, "Asin");
    NanPayload(Math.Atan, "Atan");
    NanPayload(Math.Ceiling, "Ceiling");
    NanPayload(Math.Cos, "Cos");
    NanPayload(Math.Exp, "Exp");
    NanPayload(Math.Floor, "Floor");
    NanPayload(delegate(double x) { return Math.Log(x); }, "Log");
    NanPayload(Math.Log10, "Log10");
    NanPayload(Math.Round, "Round");
    NanPayload(Math.Sin, "Sin");
    NanPayload(Math.Sqrt, "Sqrt");
    NanPayload(Math.Tan, "Tan");
    // Arithmetic operations
    NanPayload(delegate(double x, double y) { return x+y; }, "+");
    NanPayload(delegate(double x, double y) { return x-y; }, "-");
    NanPayload(delegate(double x, double y) { return x*y; }, "*");
    NanPayload(delegate(double x, double y) { return x/y; }, "/");
    NanPayload(delegate(double x, double y) { return x%y; }, "%");
    // Two-argument Math functions
    NanPayload(Math.Atan2, "Atan2");
    NanPayload(Math.IEEERemainder, "IEEERemainder");
    NanPayload(Math.Max, "Max");
    NanPayload(Math.Min, "Min");
    NanPayload(Math.Pow, "Pow");
    NanPayload(delegate(double x, double y) { return Math.Log(x,y); }, "Log");
  }

  // ------------------------------------------------------------
  // Manipulating NaNs 

  public static int ErrorCode(double d)
  {
    return (int)System.BitConverter.DoubleToInt64Bits(d);
  }

  public static double MakeNaN(int errorNumber)
  {
    long nanbits = System.BitConverter.DoubleToInt64Bits(Double.NaN);
    return System.BitConverter.Int64BitsToDouble(nanbits | errorNumber);
  }

  // ------------------------------------------------------------
  // Linear regression computed two ways

  public class Point {
    public double x, y;

    public Point(double x, double y) {
      this.x = x;
      this.y = y;
    }

    public Point move(double dx, double dy) {
      return new Point(x + dx, y + dy);
    }

    public static Point[] moveAll(Point[] ps, double dx, double dy) {
      Point[] res = new Point[ps.Length];
      for (int i=0; i<ps.Length; i++)
	res[i] = ps[i].move(dx, dy);
      return res;
    }
  }

  public static void linear1(Point[] ps) {
    int n = ps.Length;
    double SX = 0.0, SY = 0.0, SSX = 0.0, SXY = 0.0;
    for (int i=0; i<n; i++) {
      Point p = ps[i];
      SX += p.x;
      SY += p.y;
      SXY += p.x * p.y;
      SSX += p.x * p.x;
    }
    double beta = (SXY - SX*SY/n)/(SSX - SX*SX/n);
    double alpha = SY/n - SX/n * beta;
    Console.WriteLine("LINREG1: y = {0} + {0} * x", alpha, beta);
  }

  public static void linear2(Point[] ps) {
    int n = ps.Length;
    double SX = 0.0, SY = 0.0;
    for (int i=0; i<n; i++) {
      Point p = ps[i];
      SX += p.x;
      SY += p.y;
    }
    double EX = SX/n, EY = SY/n;
    double SDXDY = 0.0, SSDX = 0.0;
    for (int i=0; i<n; i++) {
      Point p = ps[i];
      double dx = p.x - EX, dy = p.y - EY;
      SDXDY += dx * dy;
      SSDX += dx * dx;
    }
    double beta = SDXDY/SSDX;
    double alpha = SY/n - SX/n * beta;
    Console.WriteLine("LINREG2: y = {0} + {0} * x", alpha, beta);
  }

  // ------------------------------------------------------------
  // From Goldberg 1991, due to Kahan; relative error only  2 eps!
  // Optimizer should not reduce C to (T-S)-Y to ((S+Y)-S)-Y to 0.0.

  public static double sumKahan(double[] xs) {
    int n = xs.Length;
    if (n == 0)
      return 0.0;
    else {
      double S = xs[0], C = 0;
      for (int i=1; i<n; i++) {
	double Y = xs[i] - C, T = S + Y;
	C = (T - S) - Y;
	S = T;
      }
      return S;
    }
  }

  // Naive summation of array of doubles:

  public static double sum(double[] xs) {
    double S = 0.0;
    for (int i=0; i<xs.Length; i++) 
      S += xs[i];
    return S;
  }

  // ------------------------------------------------------------
  // Example from Goldberg 1991 of infinities; 
  // compute x/(x^2+1) as 1/(x + 1/x) to avoid overflows; 
  // this works for x==0 too:

  public static double G(double x) {
    return 1/(x + 1/x);
  }

  // ------------------------------------------------------------
  // Example from William Kahan 1981 "Why do we need..."

  public static double R(double z) {
    return 7-3/(z-2-1/(z-7+10/(z-2-2/(z-3))));
  }

  // ------------------------------------------------------------
  // Tables of arithmetic operations on doubles, in LaTeX

  public static void AllDoubleOperations() {
    Console.WriteLine("\\subsection{Arithmetic operators}");
    Console.WriteLine(@"These all agree with the corresponding operators in Java.\\[1ex]");
    Console.WriteLine("\\tt");
    displayAll2("+", (x, y) => fix(x+y));
    displayAll2("-", (x, y) => fix(x-y));
    displayAll2("*", (x, y) => fix(x*y));
    displayAll2("/", (x, y) => fix(x/y));
    displayAll2("\\%", (x, y) => fix(x%y));
    Console.WriteLine("\\subsection{Comparison operators}");
    displayAll2("==", (x, y) => (x==y) + "");
    displayAll2("!=", (x, y) => (x!=y) + "");
    displayAll2("<", (x, y) => (x<y) + "");
    displayAll2("<=", (x, y) => (x<=y) + "");
    displayAll2(">", (x, y) =>  (x>y) + "");
    displayAll2(">=", (x, y) => (x>=y) + "");
    Console.WriteLine("\\subsection{Two-argument mathematical functions}");
    Console.WriteLine(@"\noindent\textrm{These agree with the corresponding
 functions in Java, except for \texttt{Atan2(+/-Inf,+/-Inf)}
and \texttt{Pow(NaN, +/-0.0)}.}\\[1ex]");
    displayAll2("Math.Atan2", 
		(x, y) => fix("{0:F3}", Math.Atan2(x,y)));
    displayAll2("Math.IEEERemainder", (x, y) => fix(Math.IEEERemainder(x,y)));
    displayAll2("Math.Max", (x, y) => fix(Math.Max(x,y)));
    displayAll2("Math.Min", (x, y) => fix(Math.Min(x,y)));
    displayAll2("Math.Pow", (x, y) => fix(Math.Pow(x,y)));
    Console.WriteLine("\\subsection{One-argument mathematical functions}");
    Console.WriteLine(@"\noindent\textrm{These all agree with the 
corresponding functions in Java.}\\[1ex]");
    displayAll1(
      new Function1("Math.Abs", Math.Abs),
      new Function1("Math.Acos", Math.Acos),
      new Function1("Math.Asin", Math.Asin),
      new Function1("Math.Atan", Math.Atan),
      new Function1("Math.Ceiling", Math.Ceiling),
      new Function1("Math.Cos", Math.Cos),
      new Function1("Math.Exp", Math.Exp),
      new Function1("Math.Floor", Math.Floor),
      new Function1("Math.Log", Math.Log),
      new Function1("Math.Log10", Math.Log10),
      new Function1("Math.Round", Math.Round),
      new Function1("Math.Sin", Math.Sin),
      new Function1("Math.Sqrt", Math.Sqrt),
      new Function1("Math.Tan", Math.Tan)
    );
  }

  class Function1 {
    public readonly String name;
    public readonly Func<double,double> fun;

    public Function1(String name, Func<double,double> fun) {
      this.name = name;
      this.fun = fun;
    }
    public String invoke(double x) {
      return fix("{0:F3}", fun(x));
    }
  }

  private readonly static double[] values 
  = { Double.NegativeInfinity, -2.0, -0.0, 
      0.0, +2.0, Double.PositiveInfinity, Double.NaN };
  
  private static void displayAll2(String name, Func<double,double,String> function) {
    values[2] = -0.0;  // Necessary because of bug in gmcs 2.2
    Console.WriteLine("\\noindent");
    Console.WriteLine("\\begin{tabular}{r|rrrrrrr}");
    Console.Write("\\multicolumn{1}{c|}{" + name + "}");
    for (int i=0; i<values.Length; i++) 
      Console.Write(" & " + fix(values[i]));
    Console.WriteLine("\\\\\\hline");
    for (int i=0; i<values.Length; i++) { 
      Console.Write(fix(values[i]));
      for (int j=0; j<values.Length; j++) { 
	double x = values[i], y = values[j];
	Console.Write(" & " + function(x, y));
      }
      Console.WriteLine("\\\\");
    }
    Console.WriteLine("\\end{tabular}\n\\vspace{0.5cm}\n");
  }

  private static void displayAll1(params Function1[] functions) {
    values[2] = -0.0;  // Necessary because of bug in gmcs 2.2
    Console.WriteLine("\\noindent");
    Console.WriteLine("\\begin{tabular}{l|rrrrrrr}");
    for (int i=0; i<values.Length; i++) 
      Console.Write(" & " + fix("{0:F3}", values[i]));
    Console.WriteLine("\\\\\\hline");
    for (int i=0; i<functions.Length; i++) { 
      Function1 function = functions[i];
      Console.Write(function.name);
      for (int j=0; j<values.Length; j++) { 
	double x = values[j];
	Console.Write(" & " + fix(function.invoke(x)));
      }
      Console.WriteLine("\\\\");
    }
    Console.WriteLine("\\end{tabular}\n\\vspace{0.5cm}\n");
  }

  private static String fix(String s) {
    if (s == "Infinity")
      return "+Inf";
    else if (s == "-Infinity")
      return "-Inf";
    else
      return s;
  }

  private static String fix(double r) {
    if (r == 0.0 && 1/r < 0.0)
      return fix("-" + r);
    else
      return fix("" + r);
  }

  private static String fix(String fmt, double r) {
    if (r == 0 && 1/r < 0)
      return fix("-" + String.Format(fmt, r));
    else
      return fix(String.Format(fmt, r));
  }

  // ------------------------------------------------------------
  // Exp function change November 2015

  private static void ExpFunction() {
    Console.WriteLine("\nChange to Math.Exp function, November 2015:"); 
    // Mail from Troels Damgaard 16 Nov 2015
    double c = 0.22211390926973425; 
    double x = -1.5045649238845;
    double calc = Math.Exp(x);
    ShowBoth(c, calc);

    // Other example, also Nov 2015
    var intensity = Math.Log(1.0 + 0.034);
    var expiry = 45.0;
    ShowBoth(0.033434776086237419, intensity);
    ShowBoth(0.22211390927058186, Math.Exp(-intensity * expiry));
  }

  static void ShowBoth(double x, double y) {
    displayDouble(x);
    displayDouble(y);
    Console.WriteLine(x == y ? "Equal" : "Not equal");
  }

  // ------------------------------------------------------------
  // Conversions to bit patterns

  public static void byteExamples() {
    Console.WriteLine("\n8-bit two's complement integers, byte:"); 
    displayByte(0);
    displayByte(1);
    displayByte(-1);
    displayByte(2);
    displayByte(-2);
    displayByte(6);
    displayByte(15);
    displayByte(127);
    displayByte(-127);
    displayByte(-128);
    displayByte(127+1);
  }

  public static void intExamples() {
    Console.WriteLine("\n32-bit two's complement integers, int:"); 
    displayInt(0);
    displayInt(1);
    displayInt(-1);
    displayInt(2);
    displayInt(-2);
    displayInt(6);
    displayInt(15);
    displayInt(Int32.MaxValue);
    displayInt(-Int32.MaxValue);
    displayInt(Int32.MinValue);
    unchecked {
      displayInt(Int32.MaxValue+1);
    }
  }

  public static void floatExamples() {
    Console.WriteLine("\nIEEE 754 32-bit, float:"); 
    displayFloat(0.0f);
    displayFloat(-0.0f);
    displayFloat(1.0f);
    displayFloat(0.5f);
    displayFloat(-118.625f);
    displayFloat(Single.PositiveInfinity);
    displayFloat(Single.NegativeInfinity);
    displayFloat(Single.NaN);
    displayFloat(0.1f);
    displayFloat(0.1f+0.1f+0.1f+0.1f+0.1f+0.1f+0.1f+0.1f+0.1f+0.1f);
  }

  public static void doubleExamples() {
    Console.WriteLine("\nIEEE 754 64-bit, double:"); 
    displayDouble(0.0);
    displayDouble(-0.0);
    displayDouble(1.0);
    displayDouble(0.5);
    displayDouble(-118.625);
    displayDouble(Double.PositiveInfinity);
    displayDouble(Double.NegativeInfinity);
    displayDouble(Double.NaN);
    displayDouble(Math.Log(0.0));
    displayDouble(Math.Log(-1));
    displayDouble(Math.Sqrt(-1));
    displayDouble(Math.Asin(1.5));
    displayDouble(0.1);
    displayDouble(0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1);
  }

  public static void displayByte(int n) {
    byte b = (byte)n;
    Console.WriteLine(toBits(b) + " = " + b);
  }

  public static void displayInt(int n) {
    Console.WriteLine(toBits(n) + " = " + n);
  }

  public static void displayFloat(float f) {
    Console.WriteLine(toBits(f) + " = " + f);
  }

  public static void displayDouble(double d) {
    Console.WriteLine("{0} = {1:R}", toBits(d), d);
  }

  public static String toBits(float f) {
    byte[] b = BitConverter.GetBytes(f);
    long bits;
    if (BitConverter.IsLittleEndian) // eg. x86
      bits = (b[3] << 24) | (b[2] << 16) | (b[1] << 8) | (b[0] << 0);
    else
      bits = (b[0] << 24) | (b[1] << 16) | (b[2] << 8) | (b[3] << 0);      
    StringBuilder sb = insertBlanks(toBits(bits, 32), 9, 1);
    return sb.ToString();
  }

  public static String toBits(double d) {
    long bits = BitConverter.DoubleToInt64Bits(d);
    StringBuilder sb = insertBlanks(toBits(bits, 64), 12, 1);
    return sb.ToString();
  }
  
  public static String toBits(byte b) {
    return toBits(b, 8);
  }
  
  public static String toBits(int n) {
    StringBuilder sb = insertBlanks(toBits(n, 32), 24, 16, 8);
    return sb.ToString();
  }
  
  public static String toBits(long n) {
    StringBuilder sb = insertBlanks(toBits(n, 64), 
  				    56, 48, 40, 32, 24, 16, 8);
    return sb.ToString();
  }
  
  private static StringBuilder insertBlanks(String s, params int[] pos) {
    StringBuilder sb = new StringBuilder(s);
    for (int i=0; i<pos.Length; i++)
      sb.Insert(pos[i], ' ');
    return sb;
  }
  
  public static String toBits(long n, int size) {
    char[] cs = new char[size];
    for (int i=1; i<=size; i++) {
      cs[size-i] = (n & 1) != 0 ? '1' : '0';
      n >>= 1;
    }
    return new String(cs);
  }
}