You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 
 
 
 

136 lines
3.5 KiB

#!/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 scritps
// 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("}")