#!/usr/bin/python3

import math

# Save the points in separate lists 
points_t1 = []  # A list of tuples (x, r)
points_t2 = []  # A list of tuples (x, r)
points_t3 = []  # A list of tuples (x, r)

######################################################
# Generate
######################################################

adj = 0.50  # Offset ratio for step of 1, estimate

def loop_even(upper, lower, step, ary):
  result = []
  for i in range((int)(upper * step), (int)(lower * step), -1):
    x = i / step
    z = math.exp(x + (adj / step))
    r = math.log((z+1)/(z-1))
    ary.append( (x, r) )

#####
# Create initial list which will be sorted desending

# Tbl-1 9:5 in linear steps of 0.5
loop_even(9.999, 4.999, 2, points_t1)

# Tbl-2 5:1 in linear steps of 1/16
loop_even(4.999, 0.999, 16, points_t2)

# Tbl-3 log steps:  1/(2^0.5), 1/(2^1), 1/(2^1.5), ...
for i in range(1, 28):
    e = (2 ** (i/2))
    x = 1.0 / e
    z = math.exp(x+ (.19 / e))
    r = math.log((z+1)/(z-1))
    points_t3.append( (x, r) )


######################################################
# Output
######################################################

print("""
// phi0.c
//
// An approximation of the function
//
// This file is generated by the gen_phi0 scripts
// Any changes should be made to that file, not this one

#include <stdint.h>

#define SI16(f) ((int32_t)(f * (1<<16)))

float phi0( float xf ) {

    int32_t x = SI16(xf);

    if (x >= SI16(10.0f)) return(0.0f);
""", end="")


##################################
# Tbl-1 9:5 in linear steps of 0.5  (9.5, 9.0, 8.5, 8.0, .., 5.0)
print("""    else {
      if (x >= SI16(5.0f)) {
        int i = 19 - (x >> 15);
        switch (i) {
""", end="")
for i in range(len(points_t1)):
    #assert(points_t1[i][0] == ((18 - i) / 2))
    print("{}case {}: return({:.9f}f); // ({:.1f})".format(
        (10*" "), i, points_t1[i][1], points_t1[i][0]))
print("{}}}".format(( 8*" ")))
print("{}}}".format(( 6*" ")))


##################################
# Tbl-2 5:1 in linear steps of 1/16  (4-15/16, 4-7/7, ... 1)
print("""      else {
        if (x >= SI16(1.0f)) {
          int i = 79 - (x >> 12);
          switch (i) {
""", end="")
for i in range(len(points_t2)):
    #assert(points_t2[i][0] == ((18 - i) / 2))
    print("{}case {}: return({:.9f}f); // ({:.4f})".format(
        (12*" "), i, points_t2[i][1], points_t2[i][0]))
print("{}}}".format((10*" ")))
print("{}}}".format(( 8*" ")))


##################################
# Tbl-3 log steps:  1/(2^0.5), 1/(2^1), 1/(2^1.5), ...
# Output as a balanced search

def prnt_cmp(x, ind):
    print("{}if (x > SI16({:.6f}f)) {{".format((" "*ind), x))

def prnt_rtn(r, ind):
    print("{}return({:.9f}f);".format((" "*ind), r))

def one_level(pts, ind, dft):
    #print("# One_Level({})".format(pts))
    mid = (int)(len(pts)/2)
    lft = pts[:mid]
    rgt = pts[mid+1:]
    x = pts[mid][0]
    r = pts[mid][1]
    #print("## {}, {}, {}".format(lft, x, rgt))
    prnt_cmp(x, ind)
    if (len(lft)): 
        one_level(lft, ind+2, r)
    else:
        prnt_rtn(r, (ind+2))
    print("{}}} else {{".format((" "*(ind))))
    if (len(rgt)): 
        one_level(rgt, ind+2, dft)
    else:
        prnt_rtn(dft, (ind+2))
    print("{}}}".format((" "*(ind))))
    

# Start recursive process
print("{}else {{".format((8*" ")))
indent = 10
one_level(points_t3, indent, 10)
print("{}}}".format((8*" ")))  # End of tb1_1
print("{}}}".format((6*" ")))  # End of tb1_2
print("{}}}".format((4*" ")))  # End of tb1_3
print("{}return(10.0f);".format((4*" ")))
print("}")