heron_k.h


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133

#ifndef HERON_K_H
#define HERON_K_H 1

#include <stdio.h>
/* Math functions usually need "-lm" as compiler argument: */
#include <math.h>
//#include <assert.h>
#include <float.h>
#include "c23-fallback.h"


/**
 ** @file
 ** @brief Implement a Heron process to approximate powers of `1/k`
 **/


/**
 ** @def FLT_RDXRDX
 ** @brief the radix base 2 of @c FLT_RADIX
 **
 ** This is needed internally for some of the code below.
 **/
#if FLT_RADIX == 2
# define FLT_RDXRDX 1
#elif FLT_RADIX == 4
# define FLT_RDXRDX 2
#elif FLT_RADIX == 8
# define FLT_RDXRDX 3
#elif FLT_RADIX == 16
# define FLT_RDXRDX 4
#else
# error "encoutered platform with unusual FLT_RADIX, please report"
/**
 ** @def FLT_RDXRDX
 ** @brief the radix base 2 of @c FLT_RADIX
 **
 ** This is needed internally for some of the code below.
 **/
# define FLT_RDXRDX something
#endif

/**
 ** @brief raise @a to the power of @a k
 **
 ** @warning manual implementation as an exercise, never use in
 ** production code
 **
 ** @pre @a k must be a power of 2.0
 **/
double expk2(double a, unsigned k) [[__unsequenced__]];

/**
 ** @brief raise @a to the power of @a k
 **
 ** @warning manual implementation as an exercise, never use in
 ** production code
 **
 ** @pre @a a must be strictly greater than 0.
 **/
double expk_rec(double a, unsigned k) [[__unsequenced__]];

/**
 ** @brief convert floating-point number to fractional and integral components
 **
 ** @warning manual implementation as an exercise, never use in
 ** production code
 **
 ** This implementation just repeatedly multiplies the input value by
 ** 2.0 or 0.5.
 **/
double frexp_np(double x, signed exp[static 1]) [[__unsequenced__]];

/**
 ** @brief compute a good estimate for the Heron process for k=-1
 ** @return x is a power of 2 such that 0.5 <= x*a < 1.0
 ** @pre Supposes that a > 0
 **/
double heron1_estimate(double a) [[__unsequenced__]];

/**
 ** @brief compute a good estimate for the Heron process for k=-1
 ** @return x is a power of 2 such that 0.5 <= x*a < 1.0
 ** @pre Supposes that 0 < a < 0.5
 **/
double heron1_estimate_dir_05(double a) [[__unsequenced__]];

/**
 ** @brief compute a good estimate for the Heron process for k=-1
 ** @return x is a power of 2 such that 0.5 <= x*a < 1.0
 ** @pre Supposes that a > 1.0
 **/
double heron1_estimate_dir_10(double a) [[__unsequenced__]];

/**
 ** @brief compute a good estimate for the Heron process for k=-1
 ** @return x is a power of 2 such that 0.5 <= x*a < 1.0
 ** @pre Supposes that a > 0.0
 **/
double heron1_estimate_dir(double a) [[__unsequenced__]];

/**
 ** @brief use the Heron process to approximate 1.0/a
 ** @pre Supposes that a > 0.0
 **/
double heron1(double a) [[__unsequenced__]];

/**
 ** @brief raise @a to the power of @a k
 **
 ** @warning manual implementation as an exercise, never use in
 ** production code
 **/
double expk(double a, signed k) [[__unsequenced__]];

/**
 ** @brief use the Heron process to approximate @a a to the
 ** power of `1/k`
 **
 ** Or in other words this computes the @f$k^{th}@f$ root of @a a.
 ** As a special feature, if @a k is `-1` it computes the
 ** multiplicative inverse of @a a.
 **
 ** @param a must be greater than `0.0`
 ** @param k should not be `0` and otherwise be between
 ** `DBL_MIN_EXP*FLT_RDXRDX` and
 ** `DBL_MAX_EXP*FLT_RDXRDX`.
 **
 ** @see FLT_RDXRDX
 **/
double heron(double a, signed k) [[__unsequenced__]];

#endif