<commit_before>#include "logicle.h"
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <float.h>
#include <math.h>

struct logicle_params
{
    double T, W, M, A;

    double a, b, c, d, f;
    double w, x0, x1, x2;

    double xTaylor;
    double *taylor;

    double *lookup;
    int bins;
};

const int TAYLOR_LENGTH = 16;

double solve (double b, double w) {
	// w == 0 means its really arcsinh
	if (w == 0)
		return b;

	// precision is the same as that of b
	double tolerance = 2 * b * DBL_EPSILON;

	// based on RTSAFE from Numerical Recipes 1st Edition
	// bracket the root
	double d_lo = 0;
	double d_hi = b;

	// bisection first step
	double d = (d_lo + d_hi) / 2;
	double last_delta = d_hi - d_lo;
	double delta;

	// evaluate the f(w,b) = 2 * (ln(d) - ln(b)) + w * (b + d)
	// and its derivative
	double f_b = -2 * log(b) + w * b;
	double f = 2 * log(d) + w * d + f_b;
	double last_f = NAN;

	for (int i = 1; i < 40; ++i)
	{
		// compute the derivative
		double df = 2 / d + w;

		// if Newton's method would step outside the bracket
		// or if it isn't converging quickly enough
		if (((d - d_hi) * df - f) * ((d - d_lo) * df - f) >= 0
			|| abs(1.9 * f) > abs(last_delta * df))
		{
			// take a bisection step
			delta = (d_hi - d_lo) / 2;
			d = d_lo + delta;
			if (d == d_lo)
				return d; // nothing changed, we're done
		}
		else
		{
			// otherwise take a Newton's method step
			delta = f / df;
			double t = d;
			d -= delta;
			if (d == t)
				return d; // nothing changed, we're done
		}
		// if we've reached the desired precision we're done
		if (abs(delta) < tolerance)
			return d;
		last_delta = delta;

		// recompute the function
		f = 2 * log(d) + w * d + f_b;
		if (f == 0 || f == last_f)
			return d; // found the root or are not going to get any closer
		last_f = f;

		// update the bracketing interval
		if (f < 0)
			d_lo = d;
		else
			d_hi = d;
	}

	// TODO: do something if we get here, did not converge, max iter exceeded
}

double seriesBiexponential (struct logicle_params p, double scale) {
	// Taylor series is around x1
	double x = scale - p.x1;
	// note that taylor[1] should be identically zero according
	// to the Logicle condition so skip it here
	double sum = p.taylor[TAYLOR_LENGTH - 1] * x;
	for (int i = TAYLOR_LENGTH - 2; i >= 2; --i)
		sum = (sum + p.taylor[i]) * x;
	return (sum * x + p.taylor[0]) * x;
}

double scale (struct logicle_params p, double value) {
	// handle true zero separately
	if (value == 0)
		return p.x1;

	// reflect negative values
	bool negative = value < 0;
	if (negative)
		value = -value;

	// initial guess at solution
	double x;
	if (value < p.f)
		// use linear approximation in the quasi linear region
		x = p.x1 + value / p.taylor[0];
	else
		// otherwise use ordinary logarithm
		x = log(value / p.a) / p.b;

	// try for double precision unless in extended range
	double tolerance = 3 * DBL_EPSILON;
	if (x > 1)
		tolerance = 3 * x * DBL_EPSILON;

	for (int i = 0; i < 40; ++i)
	{
		// compute the function and its first two derivatives
		double ae2bx = p.a * exp(p.b * x);
		double ce2mdx = p.c / exp(p.d * x);
		double y;
		if (x < p.xTaylor)
			// near zero use the Taylor series
			y = seriesBiexponential(p, x) - value;
		else
			// this formulation has better round-off behavior
			y = (ae2bx + p.f) - (ce2mdx + value);
		double abe2bx = p.b * ae2bx;
		double cde2mdx = p.d * ce2mdx;
		double dy = abe2bx + cde2mdx;
		double ddy = p.b * abe2bx - p.d * cde2mdx;

		// this is Halley's method with cubic convergence
		double delta = y / (dy * (1 - y * ddy / (2 * dy * dy)));
		x -= delta;

		// if we've reached the desired precision we're done
		if (abs(delta) < tolerance) {
			// handle negative arguments
			if (negative)
				return 2 * p.x1 - x;
			else
				return x;
		}
	}

	// TODO: do something if we get here, scale did not converge
};


void logicle_scale(double T, double W, double M, double A, double* x, int n) {
	// allocate the parameter structure
	struct logicle_params p;
	p.taylor = 0;

    // TODO: move these checks to Python
//	if (T <= 0)
//		throw IllegalParameter("T is not positive");
//	if (W <= 0)
//		throw IllegalParameter("W is not positive");
//	if (M <= 0)
//		throw IllegalParameter("M is not positive");
//	if (2 * W > M)
//		throw IllegalParameter("W is too large");
//	if (-A > W || A + W > M - W)<commit_msg>BUGFIX: use fabs in logicle C extension, fixes test_logicle_range test<commit_after>#include "logicle.h"
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <float.h>
#include <math.h>

struct logicle_params
{
    double T, W, M, A;

    double a, b, c, d, f;
    double w, x0, x1, x2;

    double xTaylor;
    double *taylor;

    double *lookup;
    int bins;
};

const int TAYLOR_LENGTH = 16;

double solve (double b, double w) {

	double delta;
	double df, t;

	// w == 0 means its really arcsinh
	if (w == 0)
		return b;

	// precision is the same as that of b
	double tolerance = 2 * b * DBL_EPSILON;

	// based on RTSAFE from Numerical Recipes 1st Edition
	// bracket the root
	double d_lo = 0;
	double d_hi = b;

	// bisection first step
	double d = (d_lo + d_hi) / 2;
	double last_delta = d_hi - d_lo;

	// evaluate the f(w,b) = 2 * (ln(d) - ln(b)) + w * (b + d)
	// and its derivative
	double f_b = -2 * log(b) + w * b;
	double f = 2 * log(d) + w * d + f_b;
	double last_f = NAN;

	for (int i = 1; i < 40; ++i)
	{
		// compute the derivative
		df = 2 / d + w;

		// if Newton's method would step outside the bracket
		// or if it isn't converging quickly enough
		if (((d - d_hi) * df - f) * ((d - d_lo) * df - f) >= 0
			|| fabs(1.9 * f) > fabs(last_delta * df))
		{
			// take a bisection step
			delta = (d_hi - d_lo) / 2;
			d = d_lo + delta;
			if (d == d_lo)
				return d; // nothing changed, we're done
		}
		else
		{
			// otherwise take a Newton's method step
			delta = f / df;
			t = d;
			d -= delta;
			if (d == t)
				return d; // nothing changed, we're done
		}
		// if we've reached the desired precision we're done
		if (fabs(delta) < tolerance)
			return d;
		last_delta = delta;

		// recompute the function
		f = 2 * log(d) + w * d + f_b;
		if (f == 0 || f == last_f)
			return d; // found the root or are not going to get any closer
		last_f = f;

		// update the bracketing interval
		if (f < 0)
			d_lo = d;
		else
			d_hi = d;
	}
}

double seriesBiexponential (struct logicle_params p, double scale) {
	// Taylor series is around x1
	double x = scale - p.x1;
	// note that taylor[1] should be identically zero according
	// to the Logicle condition so skip it here
	double sum = p.taylor[TAYLOR_LENGTH - 1] * x;
	for (int i = TAYLOR_LENGTH - 2; i >= 2; --i)
		sum = (sum + p.taylor[i]) * x;
	return (sum * x + p.taylor[0]) * x;
}

double scale (struct logicle_params p, double value) {
	// handle true zero separately
	if (value == 0)
		return p.x1;

	// reflect negative values
	bool negative = value < 0;
	if (negative)
		value = -value;

	// initial guess at solution
	double x;
	if (value < p.f)
		// use linear approximation in the quasi linear region
		x = p.x1 + value / p.taylor[0];
	else
		// otherwise use ordinary logarithm
		x = log(value / p.a) / p.b;

	// try for double precision unless in extended range
	double tolerance = 3 * DBL_EPSILON;
	if (x > 1)
		tolerance = 3 * x * DBL_EPSILON;

	for (int i = 0; i < 40; ++i)
	{
		// compute the function and its first two derivatives
		double ae2bx = p.a * exp(p.b * x);
		double ce2mdx = p.c / exp(p.d * x);
		double y;
		if (x < p.xTaylor)
			// near zero use the Taylor series
			y = seriesBiexponential(p, x) - value;
		else
			// this formulation has better round-off behavior
			y = (ae2bx + p.f) - (ce2mdx + value);
		double abe2bx = p.b * ae2bx;
		double cde2mdx = p.d * ce2mdx;
		double dy = abe2bx + cde2mdx;
		double ddy = p.b * abe2bx - p.d * cde2mdx;

		// this is Halley's method with cubic convergence
		double delta = y / (dy * (1 - y * ddy / (2 * dy * dy)));
		x -= delta;

		// if we've reached the desired precision we're done
		if (fabs(delta) < tolerance) {
			// handle negative arguments
			if (negative)
				return 2 * p.x1 - x;
			else
				return x;
		}
	}

	// TODO: do something if we get here, scale did not converge
};


void logicle_scale(double T, double W, double M, double A, double* x, int n) {
	// allocate the parameter structure
	struct logicle_params p;
	p.taylor = 0;

    // TODO: move these checks to Python
//	if (T <= 0)
//		throw IllegalParameter("T is not positive");
//	if (W <= 0)
//		throw IllegalParameter("W is not positive");
//	if (M <= 0)
//		throw IllegalParameter("M is not positive");
//	if (2 * W > M)
//		throw IllegalParameter("W is too large");
//	if (-A > W || A + W > M - W)
{'fraction_non_alphanumeric': 0.13473466574776016, 'fraction_numerical': 0.014702504020215943, 'mean_word_length': 2.4689243027888446, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '951844', 'n_tokens_mistral': 3690, 'n_tokens_neox': 3287, 'n_words': 1214}