diff options
Diffstat (limited to 'source/luametatex/source/mp/mpw')
| -rw-r--r-- | source/luametatex/source/mp/mpw/mp.w | 19 | ||||
| -rw-r--r-- | source/luametatex/source/mp/mpw/mpmathdecimal.w | 288 | ||||
| -rw-r--r-- | source/luametatex/source/mp/mpw/mpmathdouble.w | 138 | ||||
| -rw-r--r-- | source/luametatex/source/mp/mpw/mpmathposit.w | 1428 |
4 files changed, 1517 insertions, 356 deletions
diff --git a/source/luametatex/source/mp/mpw/mp.w b/source/luametatex/source/mp/mpw/mp.w index f4bdb58ca..7390181ea 100644 --- a/source/luametatex/source/mp/mpw/mp.w +++ b/source/luametatex/source/mp/mpw/mp.w @@ -200,6 +200,7 @@ abstraction. # include "avl.h" # include "auxmemory.h" +# include "auxposit.h" # include <string.h> # include <setjmp.h> @@ -231,6 +232,7 @@ typedef struct MP_instance { # include "mpmathdouble.h" # include "mpmathbinary.h" # include "mpmathdecimal.h" +# include "mpmathposit.h" # include "mpstrings.h" @h @<Declarations@> @@ -293,13 +295,15 @@ typedef enum mp_number_type { mp_angle_type, mp_double_type, mp_binary_type, - mp_decimal_type + mp_decimal_type, + mp_posit_type } mp_number_type; typedef union mp_number_store { void *num; double dval; int val; + posit_t pval; } mp_number_store; typedef struct mp_number_data { @@ -568,6 +572,9 @@ MP mp_initialize (MP_options * opt) case mp_math_binary_mode: mp->math = mp_initialize_binary_math(mp); break; + case mp_math_posit_mode: + mp->math = mp_initialize_posit_math(mp); + break; default: mp->math = mp_initialize_double_math(mp); break; @@ -586,6 +593,9 @@ MP mp_initialize (MP_options * opt) case mp_math_decimal_mode: set_internal_string(mp_number_system_internal, mp_intern(mp, "decimal")); break; + case mp_math_posit_mode: + set_internal_string(mp_number_system_internal, mp_intern(mp, "posit")); + break; case mp_math_binary_mode: set_internal_string(mp_number_system_internal, mp_intern(mp, "binary")); break; @@ -1868,7 +1878,8 @@ typedef enum mp_math_mode { mp_math_scaled_mode, mp_math_double_mode, mp_math_binary_mode, - mp_math_decimal_mode + mp_math_decimal_mode, + mp_math_posit_mode } mp_math_mode; @ @<Option variables@>= @@ -22558,7 +22569,7 @@ static void mp_do_unary (MP mp, int c) case mp_cos_d_operation: /* This is rather inefficient, esp decimal, to calculate both each time. We could - pass NULL as signal to do only one. + pass NULL as signal to do only one, or just have n_sin and n_cos. */ if (mp->cur_exp.type != mp_known_type) { mp_bad_unary(mp, c); @@ -22569,7 +22580,7 @@ static void mp_do_unary (MP mp, int c) new_fraction(n_sin); new_fraction(n_cos); number_clone(arg1, cur_exp_value_number); - number_clone(arg2, unity_t); + number_clone(arg2, unity_t); /* maybe dp360 */ number_multiply_int(arg2, 360); number_modulo(arg1, arg2); convert_scaled_to_angle(arg1); diff --git a/source/luametatex/source/mp/mpw/mpmathdecimal.w b/source/luametatex/source/mp/mpw/mpmathdecimal.w index 124c77ecb..359223271 100644 --- a/source/luametatex/source/mp/mpw/mpmathdecimal.w +++ b/source/luametatex/source/mp/mpw/mpmathdecimal.w @@ -227,6 +227,8 @@ mp_decimal_info mp_decimal_data = { .initialized = 0, }; +@ See mpmathdouble for documentation. @c + static void checkZero(decNumber *ret) { if (decNumberIsZero(ret) && decNumberIsNegative(ret)) { @@ -277,21 +279,26 @@ static double decNumberToDouble(decNumber *A) } } -@ Borrowed code from libdfp: +/*tex -$$ \arctan(x) = x - \frac {x^3}{3} + \frac {x^5{5} - \frac {x^7}{7} + \ldots$$ + \startformula + \arctan(x) = x - \frac {x^3}{3} + \frac {x^5{5} - \frac {x^7}{7} + \ldots + \stopformula -This power series works well, if $x$ is close to zero ($|x|<0.5$). If x is -larger, the series converges too slowly, so in order to get a smaller x, we apply -the identity + This power series works well, if $x$ is close to zero ($|x|<0.5$). If x is + larger, the series converges too slowly, so in order to get a smaller x, we apply + the identity -$$ \arctan(x) = 2 \arctan \left (\frac {\sqrt{1 + x^2}-1} {x} \right) $$ + \startformula + \arctan(x) = 2 \arctan \left (\frac {\sqrt{1 + x^2}-1} {x} \right) + \stopformula -twice. The first application gives us a new $x$ with $x < 1$. The second -application gives us a new x with $x < 0.4142136$. For that $x$, we use the power -series and multiply the result by four. + twice. The first application gives us a new $x$ with $x < 1$. The second + application gives us a new x with $x < 0.4142136$. For that $x$, we use the power + series and multiply the result by four. + +*/ -@c static void decNumberAtan(decNumber *result, decNumber *x_orig, decContext *localset) { decNumber x; @@ -373,7 +380,6 @@ static void decNumberAtan2(decNumber *result, decNumber *y, decNumber *x, decCon } } -@ @c math_data *mp_initialize_decimal_math (MP mp) { math_data *math = (math_data *) mp_memory_allocate(sizeof(math_data)); @@ -613,10 +619,6 @@ void mp_free_decimal_math (MP mp) mp_memory_free(mp->math); } -@ Creating and destruction of |mp_number| objects. Let's hope that mimalloc keeps -a pool for these. - -@ @c void mp_allocate_number (MP mp, mp_number *n, mp_number_type t) { (void) mp; @@ -625,7 +627,6 @@ void mp_allocate_number (MP mp, mp_number *n, mp_number_type t) decNumberZero(n->data.num); } -@ @c void mp_allocate_clone (MP mp, mp_number *n, mp_number_type t, mp_number *v) { (void) mp; @@ -635,7 +636,6 @@ void mp_allocate_clone (MP mp, mp_number *n, mp_number_type t, mp_number *v) decNumberCopy(n->data.num, v->data.num); } -@ @c void mp_allocate_abs (MP mp, mp_number *n, mp_number_type t, mp_number *v) { (void) mp; @@ -654,7 +654,6 @@ void mp_allocate_double (MP mp, mp_number *n, double v) decNumberFromDouble(n->data.num, v); } -@ @c void mp_free_number (MP mp, mp_number *n) { (void) mp; @@ -857,13 +856,6 @@ void mp_number_scaled_to_angle(mp_number *A) decNumberMultiply(A->data.num, A->data.num, &mp_decimal_data.angle_multiplier_decNumber, &mp_decimal_data.set); } -@* Query functions. - -@ Convert a number to a scaled value. |decNumberToInt32| is not able to make this -conversion properly, so instead we are using |decNumberToDouble| and a typecast. -Bad! - -@c int mp_number_to_scaled(mp_number *A) { int result; @@ -958,26 +950,6 @@ int mp_number_nonequalabs(mp_number *A, mp_number *B) return ! decNumberIsZero(&res); } -@ Fixed-point arithmetic is done on {\sl scaled integers} that are multiples of -$2^{-16}$. In other words, a binary point is assumed to be sixteen bit positions -from the right end of a binary computer word. - -@ One of \MP's most common operations is the calculation of -$\lfloor{a+b\over2}\rfloor$, the midpoint of two given integers |a| and~|b|. The -most decent way to do this is to write |(a+b)/2|; but on many machines it is -more efficient to calculate |(a+b)>>1|. - -Therefore the midpoint operation will always be denoted by |half(a+b)| in this -program. If \MP\ is being implemented with languages that permit binary shifting, -the |half| macro should be changed to make this operation as efficient as -possible. Since some systems have shift operators that can only be trusted to -work on positive numbers, there is also a macro |halfp| that is used only when -the quantity being halved is known to be positive or zero. - -@ Here is a procedure analogous to |print_int|. The current version is fairly -stupid, and it is not round-trip safe, but this is good enough for a beta test. - -@c char *mp_decnumber_tostring(decNumber *n) { decNumber corrected; @@ -994,7 +966,6 @@ char *mp_decimal_number_tostring (MP mp, mp_number *n) return mp_decnumber_tostring(n->data.num); } -@ @c void mp_decimal_print_number (MP mp, mp_number *n) { char *str = mp_decnumber_tostring(n->data.num); @@ -1002,51 +973,12 @@ void mp_decimal_print_number (MP mp, mp_number *n) mp_memory_free(str); } -@ Addition is not always checked to make sure that it doesn't overflow, but in -places where overflow isn't too unlikely the |slow_add| routine is used. - -@c void mp_decimal_slow_add (MP mp, mp_number *ret, mp_number *A, mp_number *B) { (void) mp; decNumberAdd(ret->data.num, A->data.num, B->data.num, &mp_decimal_data.set); } -@ The |make_fraction| routine produces the |fraction| equivalent of |p/q|, given -integers |p| and~|q|; it computes the integer -$f=\lfloor2^{28}p/q+{1\over2}\rfloor$, when $p$ and $q$ are positive. If |p| and -|q| are both of the same scaled type |t|, the \quote {type relation} -|make_fraction(t,t)=fraction| is valid; and it's also possible to use the -subroutine \quote {backwards,} using the relation |make_fraction(t,fraction)=t| -between scaled types. - -If the result would have magnitude $2^{31}$ or more, |make_fraction| sets -|arith_error:=true|. Most of \MP's internal computations have been designed to -avoid this sort of error. - -If this subroutine were programmed in assembly language on a typical machine, we -could simply compute |(@t$2^{28}$@>*p)div q|, since a double-precision product -can often be input to a fixed-point division instruction. But when we are -restricted to int-eger arithmetic it is necessary either to resort to -multiple-precision maneuvering or to use a simple but slow iteration. The -multiple-precision technique would be about three times faster than the code -adopted here, but it would be comparatively long and tricky, involving about -sixteen additional multiplications and divisions. - -This operation is part of \MP's \quote {inner loop}; indeed, it will consume nearly -10\pct! of the running time (exclusive of input and output) if the code below is -left unchanged. A machine-dependent recoding will therefore make \MP\ run faster. -The present implementation is highly portable, but slow; it avoids multiplication -and division except in the initial stage. System wizards should be careful to -replace it with a routine that is guaranteed to produce identical results in all -cases. @^system dependencies@> - -As noted below, a few more routines should also be replaced by machine-dependent -code, for efficiency. But when a procedure is not part of the \quote {inner loop,} -such changes aren't advisable; simplicity and robustness are preferable to -trickery, unless the cost is too high. @^inner loop@> - -@c void mp_decimal_make_fraction (MP mp, decNumber *ret, decNumber *p, decNumber *q) { decNumberDivide(ret, p, q, &mp_decimal_data.set); @@ -1059,16 +991,6 @@ void mp_decimal_number_make_fraction (MP mp, mp_number *ret, mp_number *p, mp_nu mp_decimal_make_fraction(mp, ret->data.num, p->data.num, q->data.num); } -@ The dual of |make_fraction| is |take_fraction|, which multiplies a given -integer~|q| by a fraction~|f|. When the operands are positive, it computes -$p=\lfloor qf/2^{28}+{1\over2}\rfloor$, a symmetric function of |q| and~|f|. - -This routine is even more \quote {inner loopy} than |make_fraction|; the present -implementation consumes almost 20\pct! of \MP's computation time during typical -jobs, so a machine-language substitute is advisable. @^inner loop@> @^system -dependencies@> - -@c void mp_decimal_take_fraction (MP mp, decNumber *ret, decNumber *p, decNumber *q) { (void) mp; @@ -1081,39 +1003,18 @@ void mp_decimal_number_take_fraction (MP mp, mp_number *ret, mp_number *p, mp_nu mp_decimal_take_fraction(mp, ret->data.num, p->data.num, q->data.num); } -@ When we want to multiply something by a |scaled| quantity, we use a scheme -analogous to |take_fraction| but with a different scaling. Given positive -operands, |take_scaled| computes the quantity $p=\lfloor -qf/2^{16}+{1\over2}\rfloor$. - -Once again it is a good idea to use a machine-language replacement if possible; -otherwise |take_scaled| will use more than 2\pct! of the running time when the -Computer Modern fonts are being generated. @^inner loop@> - -@c void mp_decimal_number_take_scaled (MP mp, mp_number *ret, mp_number *p_orig, mp_number *q_orig) { (void) mp; decNumberMultiply(ret->data.num, p_orig->data.num, q_orig->data.num, &mp_decimal_data.set); } -@ For completeness, there's also |make_scaled|, which computes a quotient as a -|scaled| number instead of as a |fraction|. In other words, the result is -$\lfloor2^{16}p/q+{1\over2}\rfloor$, if the operands are positive. \ (This -procedure is not used especially often, so it is not part of \MP's inner loop.) - -@c void mp_decimal_number_make_scaled (MP mp, mp_number *ret, mp_number *p_orig, mp_number *q_orig) { decNumberDivide(ret->data.num, p_orig->data.num, q_orig->data.num, &mp_decimal_data.set); mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set); } -@ @* Scanning numbers in the input. - -The definitions below are temporarily here - -@ @c void mp_wrapup_numeric_token (MP mp, unsigned char *start, unsigned char *stop) { decNumber result; @@ -1161,7 +1062,6 @@ void mp_wrapup_numeric_token (MP mp, unsigned char *start, unsigned char *stop) set_cur_cmd((mp_variable_type) mp_numeric_command); } -@ @c static void find_exponent (MP mp) { if (mp->buffer[mp->cur_input.loc_field] == 'e' @@ -1196,9 +1096,6 @@ void mp_decimal_scan_fractional_token (MP mp, int n) mp_wrapup_numeric_token(mp, start, stop); } -@ We just have to collect bytes. - -@c void mp_decimal_scan_numeric_token (MP mp, int n) { unsigned char *start = &mp->buffer[mp->cur_input.loc_field -1]; @@ -1218,32 +1115,6 @@ void mp_decimal_scan_numeric_token (MP mp, int n) mp_wrapup_numeric_token(mp, start, stop); } -@ The |scaled| quantities in \MP\ programs are generally supposed to be less than -$2^{12}$ in absolute value, so \MP\ does much of its internal arithmetic with -28~significant bits of precision. A |fraction| denotes a scaled integer whose -binary point is assumed to be 28 bit positions from the right. - -@ Here is a typical example of how the routines above can be used. It computes -the function $${1\over3\tau}f(\theta,\phi)= -{\tau^{-1}\bigl(2+\sqrt2\,(\sin\theta-{1\over16}\sin\phi) -(\sin\phi-{1\over16}\sin\theta)(\cos\theta-\cos\phi)\bigr)\over -3\,\bigl(1+{1\over2}(\sqrt5-1)\cos\theta+{1\over2}(3-\sqrt5\,)\cos\phi\bigr)},$$ -where $\tau$ is a |scaled| \quote {tension} parameter. This is \MP's magic fudge -factor for placing the first control point of a curve that starts at an angle -$\theta$ and ends at an angle $\phi$ from the straight path. (Actually, if the -stated quantity exceeds 4, \MP\ reduces it to~4.) - -The trigonometric quantity to be multiplied by $\sqrt2$ is less than $\sqrt2$. -(It's a sum of eight terms whose absolute values can be bounded using relations -such as $\sin\theta\cos\theta|1\over2|$.) Thus the numerator is positive; and -since the tension $\tau$ is constrained to be at least $3\over4$, the numerator -is less than $16\over3$. The denominator is nonnegative and at most~6. - -The angles $\theta$ and $\phi$ are given implicitly in terms of |fraction| -arguments |st|, |ct|, |sf|, and |cf|, representing $\sin\theta$, $\cos\theta$, -$\sin\phi$, and $\cos\phi$, respectively. - -@c void mp_decimal_velocity (MP mp, mp_number *ret, mp_number *st, mp_number *ct, mp_number *sf, mp_number *cf, mp_number *t) { decNumber acc, num, denom; /* registers for intermediate calculations */ @@ -1300,12 +1171,6 @@ void mp_decimal_velocity (MP mp, mp_number *ret, mp_number *st, mp_number *ct, m mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set); } -@ The following somewhat different subroutine tests rigorously if $ab$ is greater -than, equal to, or less than~$cd$, given integers $(a,b,c,d)$. In most cases a -quick decision is reached. The result is $+1$, 0, or~$-1$ in the three respective -cases. - -@c int mp_ab_vs_cd (mp_number *a_orig, mp_number *b_orig, mp_number *c_orig, mp_number *d_orig) { decNumber a, b, c, d; @@ -1325,33 +1190,6 @@ int mp_ab_vs_cd (mp_number *a_orig, mp_number *b_orig, mp_number *c_orig, mp_num } } -@ Now here's a subroutine that's handy for all sorts of path computations: Given -a quadratic polynomial $B(a,b,c;t)$, the |crossing_point| function returns the -unique |fraction| value |t| between 0 and~1 at which $B(a,b,c;t)$ changes from -positive to negative, or returns |t=fraction_one+1| if no such value exists. If -|a<0| (so that $B(a,b,c;t)$ is already negative at |t=0|), |crossing_point| -returns the value zero. - -The general bisection method is quite simple when $n=2$, hence |crossing_point| -does not take much time. At each stage in the recursion we have a subinterval -defined by |l| and~|j| such that $B(a,b,c;2^{-l}(j+t))=B(x_0,x_1,x_2;t)$, and we -want to \quote {zero in} on the subinterval where $x_0\G0$ and $\min(x_1,x_2)<0$. - -It is convenient for purposes of calculation to combine the values of |l| and~|j| -in a single variable $d=2^l+j$, because the operation of bisection then -corresponds simply to doubling $d$ and possibly adding~1. Furthermore it proves -to be convenient to modify our previous conventions for bisection slightly, -maintaining the variables $X_0=2^lx_0$, $X_1=2^l(x_0-x_1)$, and -$X_2=2^l(x_1-x_2)$. With these variables the conditions $x_0\ge0$ and -$\min(x_1,x_2)<0$ are equivalent to $\max(X_1,X_1+X_2)>X_0\ge0$. - -The following code maintains the invariant relations -$0\L|x0|<\max(|x1|,|x1|+|x2|)$, $\vert|x1|\vert<2^{30}$, $\vert|x2|\vert<2^{30}$; -it has been constructed in such a way that no arithmetic overflow will occur if -the inputs satisfy $a<2^{30}$, $\vert a-b\vert<2^{30}$, and $\vert -b-c\vert<2^{30}$. - -@c static void mp_decimal_crossing_point (MP mp, mp_number *ret, mp_number *aa, mp_number *bb, mp_number *cc) { decNumber a, b, c; @@ -1426,20 +1264,11 @@ static void mp_decimal_crossing_point (MP mp, mp_number *ret, mp_number *aa, mp_ mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set); } -@ We conclude this set of elementary routines with some simple rounding and -truncation operations. - -@ |round_unscaled| rounds a |scaled| and converts it to |int| - -@c int mp_round_unscaled(mp_number *x_orig) { return (int) lround(mp_number_to_double(x_orig)); } -@ |number_floor| floors a number - -@c void mp_number_floor(mp_number *i) { int round = mp_decimal_data.set.round; @@ -1448,20 +1277,12 @@ void mp_number_floor(mp_number *i) mp_decimal_data.set.round = round; } -@ |fraction_to_scaled| rounds a |fraction| and converts it to |scaled| - -@c void mp_decimal_fraction_to_round_scaled(mp_number *x_orig) { x_orig->type = mp_scaled_type; decNumberDivide(x_orig->data.num, x_orig->data.num, &mp_decimal_data.fraction_multiplier_decNumber, &mp_decimal_data.set); } -@* Algebraic and transcendental functions. \MP\ computes all of the necessary -special functions from scratch, without relying on |real| arithmetic or system -subroutines for sines, cosines, etc. - -@ @c void mp_decimal_square_rt (MP mp, mp_number *ret, mp_number *x_orig) { decNumber x; @@ -1487,9 +1308,6 @@ void mp_decimal_square_rt (MP mp, mp_number *ret, mp_number *x_orig) mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set); } -@ Pythagorean addition $\psqrt{a^2+b^2}$ is implemented by a quick hack - -@c void mp_decimal_pyth_add (MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig) { decNumber a, b; @@ -1509,9 +1327,6 @@ void mp_decimal_pyth_add (MP mp, mp_number *ret, mp_number *a_orig, mp_number *b mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set); } -@ Here is a similar algorithm for $\psqrt{a^2-b^2}$. Same quick hack, also. - -@c void mp_decimal_pyth_sub (MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig) { decNumber a, b; @@ -1545,19 +1360,12 @@ void mp_decimal_pyth_sub (MP mp, mp_number *ret, mp_number *a_orig, mp_number *b mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set); } -@ Power $a^b}$: - -@c void mp_decimal_power_of (MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig) { decNumberPower(ret->data.num, a_orig->data.num, b_orig->data.num, &mp_decimal_data.set); mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set); } -@ Here is the routine that calculates $2^8$ times the natural logarithm of a -|scaled| quantity; - -@c void mp_decimal_m_log (MP mp, mp_number *ret, mp_number *x_orig) { if (! decNumberIsPositive((decNumber *) x_orig->data.num)) { @@ -1583,10 +1391,6 @@ void mp_decimal_m_log (MP mp, mp_number *ret, mp_number *x_orig) mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set); } -@ Conversely, the exponential routine calculates $\exp(x/2^8)$, when |x| is -|scaled|. - -@c void mp_decimal_m_exp (MP mp, mp_number *ret, mp_number *x_orig) { decNumber temp, twofivesix; @@ -1606,10 +1410,6 @@ void mp_decimal_m_exp (MP mp, mp_number *ret, mp_number *x_orig) mp_decimal_data.limitedset.status = 0; } -@ Given integers |x| and |y|, not both zero, the |n_arg| function returns the -|angle| whose tangent points in the direction $(x,y)$. - -@c void mp_decimal_n_arg (MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_orig) { if (decNumberIsZero((decNumber *) x_orig->data.num) && decNumberIsZero((decNumber *) y_orig->data.num)) { @@ -1635,14 +1435,6 @@ void mp_decimal_n_arg (MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_or mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set); } -@ Conversely, the |n_sin_cos| routine takes an |angle| and produces the sine and -cosine of that angle. The results of this routine are stored in global integer -variables |n_sin| and |n_cos|. - -First, we need a decNumber function that calculates sines and cosines using the -Taylor series. This function is fairly optimized. - -@c static void sinecosine(decNumber *theangle, decNumber *c, decNumber *s) { int prec = mp_decimal_data.set.digits/2; @@ -1688,9 +1480,6 @@ static void sinecosine(decNumber *theangle, decNumber *c, decNumber *s) } } -@ Calculate sines and cosines. - -@c void mp_decimal_sin_cos (MP mp, mp_number *z_orig, mp_number *n_cos, mp_number *n_sin) { decNumber rad; @@ -1717,9 +1506,6 @@ void mp_decimal_sin_cos (MP mp, mp_number *z_orig, mp_number *n_cos, mp_number * mp_decnumber_check(mp, n_sin->data.num, &mp_decimal_data.set); } -@ This is the {\tt http://www-cs-faculty.stanford.edu/~uno/programs/rng.c} -with small cosmetic modifications. - @c # define KK 100 /* the long lag */ # define LL 37 /* the short lag */ @@ -1743,10 +1529,6 @@ static mp_decimal_random_info mp_decimal_random_data = { .ptr = &mp_decimal_random_data.dummy }; -/* put n new random numbers in aa */ -/* long aa[] destination */ -/* int n array length (must be at least KK) */ - static void ran_array(long aa[],int n) { int i, j; @@ -1764,14 +1546,6 @@ static void ran_array(long aa[],int n) } } -/* - the following routines are from exercise 3.6--15L after calling |ran_start|, - get new randoms by, e.g., "|x=ran_arr_next()|" - - Do this before using |ran_array|, |long seed| selector for different - streams. -*/ - static void ran_start(long seed) { int t, j; @@ -1826,8 +1600,6 @@ static void ran_start(long seed) mp_decimal_random_data.ptr = &mp_decimal_random_data.started; } -# define ran_arr_next() (*mp_decimal_random_data.ptr>=0? *mp_decimal_random_data.ptr++: ran_arr_cycle()) - static long ran_arr_cycle(void) { if (mp_decimal_random_data.ptr == &mp_decimal_random_data.dummy) { @@ -1840,9 +1612,6 @@ static long ran_arr_cycle(void) return mp_decimal_random_data.buf[0]; } -@ To initialize the |randoms| table, we call the following routine. - -@c void mp_init_randoms (MP mp, int seed) { int k = 1; /* more or less random integers */ @@ -1866,21 +1635,16 @@ void mp_init_randoms (MP mp, int seed) ran_start((unsigned long) seed); } -@ @c void mp_decimal_number_modulo(mp_number *a, mp_number *b) { decNumberRemainder(a->data.num, a->data.num, b->data.num, &mp_decimal_data.set); } -@ To consume a random integer for the uniform generator, the program below will -say |next_unif_random|. - -@c static void mp_next_unif_random (MP mp, mp_number *ret) { decNumber a; decNumber b; - unsigned long int op = (unsigned)ran_arr_next(); + unsigned long int op = (unsigned) (*mp_decimal_random_data.ptr>=0? *mp_decimal_random_data.ptr++: ran_arr_cycle()); (void) mp; decNumberFromInt32(&a, op); decNumberFromInt32(&b, MM); @@ -1889,9 +1653,6 @@ static void mp_next_unif_random (MP mp, mp_number *ret) mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set); } -@ To consume a random fraction, the program below will say |next_random|. - -@c static void mp_next_random (MP mp, mp_number *ret) { if (mp->j_random == 0) { @@ -1902,14 +1663,6 @@ static void mp_next_random (MP mp, mp_number *ret) mp_number_clone(ret, &(mp->randoms[mp->j_random])); } -@ To produce a uniform random number in the range |0<=u<x| or |0>=u>x| or -|0=u=x|, given a |scaled| value~|x|, we proceed as shown here. - -Note that the call of |take_fraction| will produce the values 0 and~|x| with -about half the probability that it will produce any other particular values -between 0 and~|x|, because it rounds its answers. - -@c static void mp_decimal_m_unif_rand (MP mp, mp_number *ret, mp_number *x_orig) { mp_number x, abs_x, u, y; /* |y| is trial value */ @@ -1932,11 +1685,6 @@ static void mp_decimal_m_unif_rand (MP mp, mp_number *ret, mp_number *x_orig) mp_free_number(mp, &u); } -@ Finally, a normal deviate with mean zero and unit standard deviation can -readily be obtained with the ratio method (Algorithm 3.4.1R in {\sl The Art of -Computer Programming}). - -@c static void mp_decimal_m_norm_rand (MP mp, mp_number *ret) { mp_number abs_x, u, r, la, xa; diff --git a/source/luametatex/source/mp/mpw/mpmathdouble.w b/source/luametatex/source/mp/mpw/mpmathdouble.w index 6f0f8df99..f5a91df75 100644 --- a/source/luametatex/source/mp/mpw/mpmathdouble.w +++ b/source/luametatex/source/mp/mpw/mpmathdouble.w @@ -71,7 +71,6 @@ First, here are some very important constants. @ Here are the functions that are static as they are not used elsewhere. @<Declarations@>= -static int mp_ab_vs_cd (mp_number *a, mp_number *b, mp_number *c, mp_number *d); static void mp_allocate_abs (MP mp, mp_number *n, mp_number_type t, mp_number *v); static void mp_allocate_clone (MP mp, mp_number *n, mp_number_type t, mp_number *v); static void mp_allocate_double (MP mp, mp_number *n, double v); @@ -306,7 +305,7 @@ math_data *mp_initialize_double_math(MP mp) math->md_print = mp_double_print_number; math->md_tostring = mp_double_number_tostring; math->md_modulo = mp_number_modulo; - math->md_ab_vs_cd = mp_ab_vs_cd; + math->md_ab_vs_cd = mp_double_ab_vs_cd; math->md_crossing_point = mp_double_crossing_point; math->md_scan_numeric = mp_double_scan_numeric_token; math->md_scan_fractional = mp_double_scan_fractional_token; @@ -353,7 +352,7 @@ void mp_free_double_math (MP mp) @ Creating an destroying |mp_number| objects -@ @c +@c void mp_allocate_number (MP mp, mp_number *n, mp_number_type t) { (void) mp; @@ -361,7 +360,6 @@ void mp_allocate_number (MP mp, mp_number *n, mp_number_type t) n->type = t; } -@ @c void mp_allocate_clone (MP mp, mp_number *n, mp_number_type t, mp_number *v) { (void) mp; @@ -369,7 +367,6 @@ void mp_allocate_clone (MP mp, mp_number *n, mp_number_type t, mp_number *v) n->data.dval = v->data.dval; } -@ @c void mp_allocate_abs (MP mp, mp_number *n, mp_number_type t, mp_number *v) { (void) mp; @@ -377,7 +374,6 @@ void mp_allocate_abs (MP mp, mp_number *n, mp_number_type t, mp_number *v) n->data.dval = fabs(v->data.dval); } -@ @c void mp_allocate_double (MP mp, mp_number *n, double v) { (void) mp; @@ -385,7 +381,6 @@ void mp_allocate_double (MP mp, mp_number *n, double v) n->data.dval = v; } -@ @c void mp_free_number (MP mp, mp_number *n) { (void) mp; @@ -901,46 +896,16 @@ quick decision is reached. The result is $+1$, 0, or~$-1$ in the three respectiv cases. @c -int mp_ab_vs_cd (mp_number *a_orig, mp_number *b_orig, mp_number *c_orig, mp_number *d_orig) +int mp_double_ab_vs_cd (mp_number *a_orig, mp_number *b_orig, mp_number *c_orig, mp_number *d_orig) { - return mp_double_ab_vs_cd(a_orig, b_orig, c_orig, d_orig); -} - -@ @<Reduce to the case that |a...@>= -if (a < 0) { - a = -a; - b = -b; -} -if (c < 0) { - c = -c; - d = -d; -} -if (d <= 0) { - if (b >= 0) { - if ((a == 0 || b == 0) && (c == 0 || d == 0)) { - ret->data.dval = 0; - } else { - ret->data.dval = 1; - } - goto RETURN; - } if (d == 0) { - ret->data.dval = (a == 0 ? 0 : -1); - goto RETURN; - } else - q = a; - a = c; - c = q; - q = -b; - b = -d; - d = q; - } -} else if (b <= 0) { - if (b < 0 && a > 0) { - ret->data.dval = -1; - return; - } else - ret->data.dval = (c == 0 ? 0 : -1); - goto RETURN; + double ab = a_orig->data.dval * b_orig->data.dval; + double cd = c_orig->data.dval * d_orig->data.dval; + if (ab > cd) { + return 1; + } else if (ab < cd) { + return -1; + } else { + return 0; } } @@ -1210,8 +1175,9 @@ void mp_double_n_arg (MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_ori } else { ret->type = mp_angle_type; ret->data.dval = atan2(y_orig->data.dval, x_orig->data.dval) * (180.0 / PI) * angle_multiplier; - if (ret->data.dval == -0.0) - ret->data.dval = 0.0; + if (ret->data.dval == -0.0) { + ret->data.dval = 0.0; + } } } @@ -1249,7 +1215,8 @@ void mp_double_sin_cos (MP mp, mp_number *z_orig, mp_number *n_cos, mp_number *n } @ This is the http://www-cs-faculty.stanford.edu/~uno/programs/rng.c with small -cosmetic modifications. +cosmetic modifications. The code only documented here as the other non scaled +number models use the same method. @c # define KK 100 /* the long lag */ @@ -1352,8 +1319,6 @@ static void mp_double_aux_ran_start(long seed) mp_double_random_data.ptr = &mp_double_random_data.started; } -# define mp_double_aux_ran_arr_next() (*mp_double_random_data.ptr>=0? *mp_double_random_data.ptr++: mp_double_aux_ran_arr_cycle()) - static long mp_double_aux_ran_arr_cycle(void) { if (mp_double_random_data.ptr == &mp_double_random_data.dummy) { @@ -1394,16 +1359,6 @@ void mp_init_randoms (MP mp, int seed) } @ Here |frac| contains what's beyond the |.|. @c -/* -static double modulus(double left, double right) -{ - double quota = left / right; - double tmp; - double frac = modf(quota, &tmp); - frac *= right; - return frac; -} -*/ void mp_number_modulo(mp_number *a, mp_number *b) { @@ -1417,7 +1372,7 @@ say |next_unif_random|. @c static void mp_next_unif_random (MP mp, mp_number *ret) { - unsigned long int op = (unsigned) mp_double_aux_ran_arr_next(); + unsigned long int op = (unsigned) (*mp_double_random_data.ptr>=0? *mp_double_random_data.ptr++: mp_double_aux_ran_arr_cycle()); double a = op / (MM * 1.0); (void) mp; ret->data.dval = a; @@ -1437,11 +1392,10 @@ static void mp_next_random (MP mp, mp_number *ret) } @ To produce a uniform random number in the range |0<=u<x| or |0>=u>x| or -|0=u=x|, given a |scaled| value~|x|, we proceed as shown here. - -Note that the call of |take_fraction| will produce the values 0 and~|x| with -about half the probability that it will produce any other particular values -between 0 and~|x|, because it rounds its answers. +|0=u=x|, given a |scaled| value~|x|, we proceed as shown here. Note that the +call of |take_fraction| will produce the values 0 and~|x| with about half +the probability that it will produce any other particular values between 0 +and~|x|, because it rounds its answers. @c static void mp_double_m_unif_rand (MP mp, mp_number *ret, mp_number *x_orig) @@ -1504,20 +1458,40 @@ static void mp_double_m_norm_rand (MP mp, mp_number *ret) mp_free_number(mp, &u); } -@ The following subroutine is used only in |norm_rand| and tests if $ab$ is -greater than, equal to, or less than~$cd$. The result is $+1$, 0, or~$-1$ in the -three respective cases. - -@c -int mp_double_ab_vs_cd (mp_number *a_orig, mp_number *b_orig, mp_number *c_orig, mp_number *d_orig) -{ - double ab = a_orig->data.dval * b_orig->data.dval; - double cd = c_orig->data.dval * d_orig->data.dval; - if (ab > cd) { - return 1; - } else if (ab < cd) { - return -1; - } else { - return 0; +@ @<Reduce to the case that |a...@>= +if (a < 0) { + a = -a; + b = -b; +} +if (c < 0) { + c = -c; + d = -d; +} +if (d <= 0) { + if (b >= 0) { + if ((a == 0 || b == 0) && (c == 0 || d == 0)) { + ret->data.dval = 0; + } else { + ret->data.dval = 1; + } + goto RETURN; + } if (d == 0) { + ret->data.dval = (a == 0 ? 0 : -1); + goto RETURN; + } else + q = a; + a = c; + c = q; + q = -b; + b = -d; + d = q; + } +} else if (b <= 0) { + if (b < 0 && a > 0) { + ret->data.dval = -1; + return; + } else + ret->data.dval = (c == 0 ? 0 : -1); + goto RETURN; } } diff --git a/source/luametatex/source/mp/mpw/mpmathposit.w b/source/luametatex/source/mp/mpw/mpmathposit.w new file mode 100644 index 000000000..8477ef27b --- /dev/null +++ b/source/luametatex/source/mp/mpw/mpmathposit.w @@ -0,0 +1,1428 @@ +% This file is part of MetaPost. The MetaPost program is in the public domain. + +@ Introduction. + +TODO: collect constants like decimal +TODO: share scanners and random + +@c +# include "mpconfig.h" +# include "mpmathposit.h" + +@h + +@ @c +@<Declarations@> + +@ @(mpmathposit.h@>= +# ifndef MPMATHPOSIT_H +# define MPMATHPOSIT_H 1 + +# include "mp.h" +# include "softposit.h" + +math_data *mp_initialize_posit_math (MP mp); + +# endif + +@* Math initialization. + +First, here are some very important constants. + +@d mp_fraction_multiplier 4096 +@d mp_angle_multiplier 16 +@d mp_warning_limit pow(2.0,52) + +@d odd(A) (abs(A)%2==1) + +@d two_to_the(A) (1<<(unsigned)(A)) +@d set_cur_cmd(A) mp->cur_mod_->command = (A) +@d set_cur_mod(A) mp->cur_mod_->data.n.data.pval = (A) + +@ Here are the functions that are static as they are not used elsewhere. + +@<Declarations@>= +static void mp_allocate_abs (MP mp, mp_number *n, mp_number_type t, mp_number *v); +static void mp_allocate_clone (MP mp, mp_number *n, mp_number_type t, mp_number *v); +static void mp_allocate_double (MP mp, mp_number *n, double v); +static void mp_allocate_number (MP mp, mp_number *n, mp_number_type t); +static int mp_posit_ab_vs_cd (mp_number *a, mp_number *b, mp_number *c, mp_number *d); +static void mp_posit_abs (mp_number *A); +static void mp_posit_crossing_point (MP mp, mp_number *ret, mp_number *a, mp_number *b, mp_number *c); +static void mp_posit_fraction_to_round_scaled (mp_number *x); +static void mp_posit_m_exp (MP mp, mp_number *ret, mp_number *x_orig); +static void mp_posit_m_log (MP mp, mp_number *ret, mp_number *x_orig); +static void mp_posit_m_norm_rand (MP mp, mp_number *ret); +static void mp_posit_m_unif_rand (MP mp, mp_number *ret, mp_number *x_orig); +static void mp_posit_n_arg (MP mp, mp_number *ret, mp_number *x, mp_number *y); +static void mp_posit_number_make_fraction (MP mp, mp_number *r, mp_number *p, mp_number *q); +static void mp_posit_number_make_scaled (MP mp, mp_number *r, mp_number *p, mp_number *q); +static void mp_posit_number_take_fraction (MP mp, mp_number *r, mp_number *p, mp_number *q); +static void mp_posit_number_take_scaled (MP mp, mp_number *r, mp_number *p, mp_number *q); +static void mp_posit_power_of (MP mp, mp_number *r, mp_number *a, mp_number *b); +static void mp_posit_print_number (MP mp, mp_number *n); +static void mp_posit_pyth_add (MP mp, mp_number *r, mp_number *a, mp_number *b); +static void mp_posit_pyth_sub (MP mp, mp_number *r, mp_number *a, mp_number *b); +static void mp_posit_scan_fractional_token (MP mp, int n); +static void mp_posit_scan_numeric_token (MP mp, int n); +static void mp_posit_set_precision (MP mp); +static void mp_posit_sin_cos (MP mp, mp_number *z_orig, mp_number *n_cos, mp_number *n_sin); +static void mp_posit_slow_add (MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_orig); +static void mp_posit_square_rt (MP mp, mp_number *ret, mp_number *x_orig); +static void mp_posit_velocity (MP mp, mp_number *ret, mp_number *st, mp_number *ct, mp_number *sf, mp_number *cf, mp_number *t); +static void mp_free_posit_math (MP mp); +static void mp_free_number (MP mp, mp_number *n); +static void mp_init_randoms (MP mp, int seed); +static void mp_number_abs_clone (mp_number *A, mp_number *B); +static void mp_number_add (mp_number *A, mp_number *B); +static void mp_number_add_scaled (mp_number *A, int B); /* also for negative B */ +static void mp_number_angle_to_scaled (mp_number *A); +static void mp_number_clone (mp_number *A, mp_number *B); +static void mp_number_divide_int (mp_number *A, int B); +static void mp_number_double (mp_number *A); +static int mp_number_equal (mp_number *A, mp_number *B); +static void mp_number_floor (mp_number *i); +static void mp_number_fraction_to_scaled (mp_number *A); +static int mp_number_greater (mp_number *A, mp_number *B); +static void mp_number_half (mp_number *A); +static int mp_number_less (mp_number *A, mp_number *B); +static void mp_number_modulo (mp_number *a, mp_number *b); +static void mp_number_multiply_int (mp_number *A, int B); +static void mp_number_negate (mp_number *A); +static void mp_number_negated_clone (mp_number *A, mp_number *B); +static int mp_number_nonequalabs (mp_number *A, mp_number *B); +static int mp_number_odd (mp_number *A); +static void mp_number_scaled_to_angle (mp_number *A); +static void mp_number_scaled_to_fraction (mp_number *A); +static void mp_number_subtract (mp_number *A, mp_number *B); +static void mp_number_swap (mp_number *A, mp_number *B); +static int mp_number_to_boolean (mp_number *A); +static double mp_number_to_double (mp_number *A); +static int mp_number_to_int (mp_number *A); +static int mp_number_to_scaled (mp_number *A); +static int mp_round_unscaled (mp_number *x_orig); +static void mp_set_posit_from_addition (mp_number *A, mp_number *B, mp_number *C); +static void mp_set_posit_from_boolean (mp_number *A, int B); +static void mp_set_posit_from_div (mp_number *A, mp_number *B, mp_number *C); +static void mp_set_posit_from_double (mp_number *A, double B); +static void mp_set_posit_from_int (mp_number *A, int B); +static void mp_set_posit_from_int_div (mp_number *A, mp_number *B, int C); +static void mp_set_posit_from_int_mul (mp_number *A, mp_number *B, int C); +static void mp_set_posit_from_mul (mp_number *A, mp_number *B, mp_number *C); +static void mp_set_posit_from_of_the_way (MP mp, mp_number *A, mp_number *t, mp_number *B, mp_number *C); +static void mp_set_posit_from_scaled (mp_number *A, int B); +static void mp_set_posit_from_subtraction (mp_number *A, mp_number *B, mp_number *C); +static void mp_set_posit_half_from_addition (mp_number *A, mp_number *B, mp_number *C); +static void mp_set_posit_half_from_subtraction (mp_number *A, mp_number *B, mp_number *C); +static void mp_wrapup_numeric_token (MP mp, unsigned char *start, unsigned char *stop); +static char *mp_posit_number_tostring (MP mp, mp_number *n); + +typedef struct mp_posit_info { + posit_t unity; + posit_t zero; + posit_t one; + posit_t two; + posit_t three; + posit_t four; + posit_t five; + posit_t eight; + posit_t seven; + posit_t sixteen; + posit_t half_unit; + posit_t minusone; + posit_t three_quarter_unit; + posit_t d16; + posit_t d64; + posit_t d256; + posit_t d4096; + posit_t d65536; + posit_t dp90; + posit_t dp180; + posit_t dp270; + posit_t dp360; + posit_t dm90; + posit_t dm180; + posit_t dm270; + posit_t dm360; + posit_t fraction_multiplier; + posit_t negative_fraction_multiplier; /* todo: also in decimal */ + posit_t angle_multiplier; + posit_t fraction_one; + posit_t fraction_two; + posit_t fraction_three; + posit_t fraction_four; + posit_t fraction_half; + posit_t fraction_one_and_half; + posit_t one_eighty_degrees; + posit_t negative_one_eighty_degrees; + posit_t three_sixty_degrees; + posit_t no_crossing; + posit_t one_crossing; + posit_t zero_crossing; + posit_t error_correction; + posit_t pi; + posit_t pi_divided_by_180; + posit_t epsilon; + posit_t EL_GORDO; + posit_t negative_EL_GORDO; + posit_t one_third_EL_GORDO; + posit_t coef; + posit_t coef_bound; + posit_t scaled_threshold; + posit_t fraction_threshold; + posit_t equation_threshold; + posit_t near_zero_angle; + posit_t p_over_v_threshold; + posit_t warning_limit; + posit_t sqrt_two_mul_fraction_one; + posit_t sqrt_five_minus_one_mul_fraction_one_and_half; + posit_t three_minus_sqrt_five_mul_fraction_one_and_half; + posit_t d180_divided_by_pi_mul_angle; + int initialized; +} mp_posit_info; + +mp_posit_info mp_posit_data = { + .initialized = 0, +}; + +inline static posit_t mp_posit_make_fraction (posit_t p, posit_t q) { return posit_mul(posit_div(p,q), mp_posit_data.fraction_multiplier); } +inline static posit_t mp_posit_take_fraction (posit_t p, posit_t q) { return posit_div(posit_mul(p,q), mp_posit_data.fraction_multiplier); } +inline static posit_t mp_posit_make_scaled (posit_t p, posit_t q) { return posit_div(p,q); } + +math_data *mp_initialize_posit_math(MP mp) +{ + math_data *math = (math_data *) mp_memory_allocate(sizeof(math_data)); + /* alloc */ + if (! mp_posit_data.initialized) { + mp_posit_data.initialized = 1; + mp_posit_data.unity = integer_to_posit(1); + mp_posit_data.zero = integer_to_posit(0); + mp_posit_data.one = integer_to_posit(1); + mp_posit_data.two = integer_to_posit(2); + mp_posit_data.three = integer_to_posit(3); + mp_posit_data.four = integer_to_posit(4); + mp_posit_data.five = integer_to_posit(5); + mp_posit_data.seven = integer_to_posit(7); + mp_posit_data.eight = integer_to_posit(8); + mp_posit_data.sixteen = integer_to_posit(16); + mp_posit_data.dp90 = integer_to_posit(90); + mp_posit_data.dp180 = integer_to_posit(180); + mp_posit_data.dp270 = integer_to_posit(270); + mp_posit_data.dp360 = integer_to_posit(360); + mp_posit_data.dm90 = integer_to_posit(-90); + mp_posit_data.dm180 = integer_to_posit(-180); + mp_posit_data.dm270 = integer_to_posit(-270); + mp_posit_data.dm360 = integer_to_posit(-360); + mp_posit_data.d16 = integer_to_posit(16); + mp_posit_data.d64 = integer_to_posit(64); + mp_posit_data.d256 = integer_to_posit(256); + mp_posit_data.d4096 = integer_to_posit(4096); + mp_posit_data.d65536 = integer_to_posit(65536); + mp_posit_data.minusone = posit_neg(mp_posit_data.one); + mp_posit_data.half_unit = posit_div(mp_posit_data.unity, mp_posit_data.two); + mp_posit_data.three_quarter_unit = posit_mul(mp_posit_data.three, posit_div(mp_posit_data.unity,mp_posit_data.four)); + mp_posit_data.fraction_multiplier = integer_to_posit(mp_fraction_multiplier); + mp_posit_data.negative_fraction_multiplier = posit_neg(mp_posit_data.fraction_multiplier); + mp_posit_data.angle_multiplier = integer_to_posit(mp_angle_multiplier); + mp_posit_data.fraction_one = mp_posit_data.fraction_multiplier; + mp_posit_data.fraction_two = posit_mul(mp_posit_data.fraction_multiplier, mp_posit_data.two); + mp_posit_data.fraction_three = posit_mul(mp_posit_data.fraction_multiplier, mp_posit_data.three); + mp_posit_data.fraction_four = posit_mul(mp_posit_data.fraction_multiplier, mp_posit_data.four); + mp_posit_data.fraction_half = posit_div(mp_posit_data.fraction_multiplier, mp_posit_data.two); + mp_posit_data.fraction_one_and_half = posit_add(mp_posit_data.fraction_multiplier, mp_posit_data.fraction_half); + mp_posit_data.one_eighty_degrees = posit_mul(mp_posit_data.angle_multiplier, mp_posit_data.dp180); + mp_posit_data.negative_one_eighty_degrees = posit_mul(mp_posit_data.angle_multiplier, mp_posit_data.dm180); + mp_posit_data.three_sixty_degrees = posit_mul(mp_posit_data.angle_multiplier, mp_posit_data.dp360); + mp_posit_data.no_crossing = posit_add(mp_posit_data.fraction_multiplier, mp_posit_data.one); + mp_posit_data.one_crossing = mp_posit_data.fraction_multiplier; + mp_posit_data.zero_crossing = mp_posit_data.zero; + mp_posit_data.error_correction = double_to_posit(1E-12); /* debatable */ + mp_posit_data.warning_limit = posit_pow(mp_posit_data.two, integer_to_posit(52)); /* this is a large value that can just be expressed without loss of precision */ + mp_posit_data.pi = double_to_posit(3.1415926535897932384626433832795028841971); + mp_posit_data.pi_divided_by_180 = posit_div(mp_posit_data.pi, mp_posit_data.dp180); + mp_posit_data.epsilon = posit_pow(mp_posit_data.two, integer_to_posit(-52.0)); + mp_posit_data.EL_GORDO = posit_sub(posit_div(double_to_posit(DBL_MAX),mp_posit_data.two), mp_posit_data.one); /* the largest value that \MP\ likes. */ + mp_posit_data.negative_EL_GORDO = posit_neg(mp_posit_data.EL_GORDO); + mp_posit_data.one_third_EL_GORDO = posit_div(mp_posit_data.EL_GORDO, mp_posit_data.three); + mp_posit_data.coef = posit_div(mp_posit_data.seven, mp_posit_data.three); /* |fraction| approximation to 7/3 */ + mp_posit_data.coef_bound = posit_mul(mp_posit_data.coef, mp_posit_data.fraction_multiplier); + mp_posit_data.scaled_threshold = double_to_posit(0.000122); /* a |scaled| coefficient less than this is zeroed */ + mp_posit_data.near_zero_angle = posit_mul(double_to_posit(0.0256), mp_posit_data.angle_multiplier); /* an angle of about 0.0256 */ + mp_posit_data.p_over_v_threshold = integer_to_posit(0x80000); + mp_posit_data.equation_threshold = double_to_posit(0.001); + + mp_posit_data.sqrt_two_mul_fraction_one = + posit_mul( + posit_sqrt(mp_posit_data.two), + mp_posit_data.fraction_one + ); + + mp_posit_data.sqrt_five_minus_one_mul_fraction_one_and_half = + posit_mul( + posit_mul( + mp_posit_data.three, + mp_posit_data.fraction_half + ), + posit_sub( + posit_sqrt(mp_posit_data.five), + mp_posit_data.one + ) + ); + + mp_posit_data.three_minus_sqrt_five_mul_fraction_one_and_half = + posit_mul( + posit_mul( + mp_posit_data.three, + mp_posit_data.fraction_half + ), + posit_sub( + mp_posit_data.three, + posit_sqrt(mp_posit_data.five) + ) + ); + + mp_posit_data.d180_divided_by_pi_mul_angle = + posit_mul( + posit_div( + mp_posit_data.dp180, + mp_posit_data.pi + ), + mp_posit_data.angle_multiplier + ); + + } + /* alloc */ + math->md_allocate = mp_allocate_number; + math->md_free = mp_free_number; + math->md_allocate_clone = mp_allocate_clone; + math->md_allocate_abs = mp_allocate_abs; + math->md_allocate_double = mp_allocate_double; + /* precission */ + mp_allocate_number(mp, &math->md_precision_default, mp_scaled_type); + mp_allocate_number(mp, &math->md_precision_max, mp_scaled_type); + mp_allocate_number(mp, &math->md_precision_min, mp_scaled_type); + /* here are the constants for |scaled| objects */ + mp_allocate_number(mp, &math->md_epsilon_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_inf_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_negative_inf_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_warning_limit_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_one_third_inf_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_unity_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_two_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_three_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_half_unit_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_three_quarter_unit_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_zero_t, mp_scaled_type); + /* |fractions| */ + mp_allocate_number(mp, &math->md_arc_tol_k, mp_fraction_type); + mp_allocate_number(mp, &math->md_fraction_one_t, mp_fraction_type); + mp_allocate_number(mp, &math->md_fraction_half_t, mp_fraction_type); + mp_allocate_number(mp, &math->md_fraction_three_t, mp_fraction_type); + mp_allocate_number(mp, &math->md_fraction_four_t, mp_fraction_type); + /* |angles| */ + mp_allocate_number(mp, &math->md_three_sixty_deg_t, mp_angle_type); + mp_allocate_number(mp, &math->md_one_eighty_deg_t, mp_angle_type); + mp_allocate_number(mp, &math->md_negative_one_eighty_deg_t, mp_angle_type); + /* various approximations */ + mp_allocate_number(mp, &math->md_one_k, mp_scaled_type); + mp_allocate_number(mp, &math->md_sqrt_8_e_k, mp_scaled_type); + mp_allocate_number(mp, &math->md_twelve_ln_2_k, mp_fraction_type); + mp_allocate_number(mp, &math->md_coef_bound_k, mp_fraction_type); + mp_allocate_number(mp, &math->md_coef_bound_minus_1, mp_fraction_type); + mp_allocate_number(mp, &math->md_twelvebits_3, mp_scaled_type); + mp_allocate_number(mp, &math->md_twentysixbits_sqrt2_t, mp_fraction_type); + mp_allocate_number(mp, &math->md_twentyeightbits_d_t, mp_fraction_type); + mp_allocate_number(mp, &math->md_twentysevenbits_sqrt2_d_t, mp_fraction_type); + /* thresholds */ + mp_allocate_number(mp, &math->md_fraction_threshold_t, mp_fraction_type); + mp_allocate_number(mp, &math->md_half_fraction_threshold_t, mp_fraction_type); + mp_allocate_number(mp, &math->md_scaled_threshold_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_half_scaled_threshold_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_near_zero_angle_t, mp_angle_type); + mp_allocate_number(mp, &math->md_p_over_v_threshold_t, mp_fraction_type); + mp_allocate_number(mp, &math->md_equation_threshold_t, mp_scaled_type); + /* initializations */ + math->md_precision_default.data.pval = posit_mul(mp_posit_data.d16, mp_posit_data.unity); + math->md_precision_max.data.pval = posit_mul(mp_posit_data.d16, mp_posit_data.unity); + math->md_precision_min.data.pval = posit_mul(mp_posit_data.d16, mp_posit_data.unity); + math->md_epsilon_t.data.pval = mp_posit_data.epsilon; + math->md_inf_t.data.pval = mp_posit_data.EL_GORDO; + math->md_negative_inf_t.data.pval = mp_posit_data.negative_EL_GORDO; + math->md_one_third_inf_t.data.pval = mp_posit_data.one_third_EL_GORDO; + math->md_warning_limit_t.data.pval = mp_posit_data.warning_limit; + math->md_unity_t.data.pval = mp_posit_data.unity; + math->md_two_t.data.pval = mp_posit_data.two; + math->md_three_t.data.pval = mp_posit_data.three; + math->md_half_unit_t.data.pval = mp_posit_data.half_unit; + math->md_three_quarter_unit_t.data.pval = mp_posit_data.three_quarter_unit; + math->md_arc_tol_k.data.pval = posit_div(mp_posit_data.unity, mp_posit_data.d4096); /* quit when change in arc length estimate reaches this */ + math->md_fraction_one_t.data.pval = mp_posit_data.fraction_one; + math->md_fraction_half_t.data.pval = mp_posit_data.fraction_half; + math->md_fraction_three_t.data.pval = mp_posit_data.fraction_three; + math->md_fraction_four_t.data.pval = mp_posit_data.fraction_four; + math->md_three_sixty_deg_t.data.pval = mp_posit_data.three_sixty_degrees; + math->md_one_eighty_deg_t.data.pval = mp_posit_data.one_eighty_degrees; + math->md_negative_one_eighty_deg_t.data.pval = mp_posit_data.negative_one_eighty_degrees; + math->md_one_k.data.pval = posit_div(mp_posit_data.one, mp_posit_data.d64); + math->md_sqrt_8_e_k.data.pval = double_to_posit(1.71552776992141359295); /* $2^{16}\sqrt{8/e} \approx 112428.82793$ */ + math->md_twelve_ln_2_k.data.pval = posit_mul(double_to_posit(8.31776616671934371292), mp_posit_data.d256); /* $2^{24}\cdot12\ln2 \approx139548959.6165 $ */ + math->md_twelvebits_3.data.pval = posit_div(integer_to_posit(1365), mp_posit_data.unity); /* $1365 \approx 2^{12}/3 $ */ + math->md_twentysixbits_sqrt2_t.data.pval = posit_div(integer_to_posit(94906266), mp_posit_data.d65536); /* $2^{26}\sqrt2 \approx 94906265.62 $ */ + math->md_twentyeightbits_d_t.data.pval = posit_div(integer_to_posit(35596755), mp_posit_data.d65536); /* $2^{28}d \approx 35596754.69 $ */ + math->md_twentysevenbits_sqrt2_d_t.data.pval = posit_div(integer_to_posit(25170707), mp_posit_data.d65536); /* $2^{27}\sqrt2\,d \approx 25170706.63 $ */ + math->md_coef_bound_k.data.pval = mp_posit_data.coef_bound; + math->md_coef_bound_minus_1.data.pval = posit_sub(mp_posit_data.coef_bound, posit_div(mp_posit_data.one, mp_posit_data.d65536)); + math->md_fraction_threshold_t.data.pval = double_to_posit(0.04096); /* a |fraction| coefficient less than this is zeroed */ + math->md_half_fraction_threshold_t.data.pval = posit_div(mp_posit_data.fraction_threshold, mp_posit_data.two); + math->md_scaled_threshold_t.data.pval = mp_posit_data.scaled_threshold; + math->md_half_scaled_threshold_t.data.pval = posit_div(mp_posit_data.scaled_threshold,mp_posit_data.two); + math->md_near_zero_angle_t.data.pval = mp_posit_data.near_zero_angle; + math->md_p_over_v_threshold_t.data.pval = mp_posit_data.p_over_v_threshold; + math->md_equation_threshold_t.data.pval = mp_posit_data.equation_threshold; + + /* functions */ + math->md_from_int = mp_set_posit_from_int; + math->md_from_boolean = mp_set_posit_from_boolean; + math->md_from_scaled = mp_set_posit_from_scaled; + math->md_from_double = mp_set_posit_from_double; + math->md_from_addition = mp_set_posit_from_addition; + math->md_half_from_addition = mp_set_posit_half_from_addition; + math->md_from_subtraction = mp_set_posit_from_subtraction; + math->md_half_from_subtraction = mp_set_posit_half_from_subtraction; + math->md_from_oftheway = mp_set_posit_from_of_the_way; + math->md_from_div = mp_set_posit_from_div; + math->md_from_mul = mp_set_posit_from_mul; + math->md_from_int_div = mp_set_posit_from_int_div; + math->md_from_int_mul = mp_set_posit_from_int_mul; + math->md_negate = mp_number_negate; + math->md_add = mp_number_add; + math->md_subtract = mp_number_subtract; + math->md_half = mp_number_half; + math->md_do_double = mp_number_double; + math->md_abs = mp_posit_abs; + math->md_clone = mp_number_clone; + math->md_negated_clone = mp_number_negated_clone; + math->md_abs_clone = mp_number_abs_clone; + math->md_swap = mp_number_swap; + math->md_add_scaled = mp_number_add_scaled; + math->md_multiply_int = mp_number_multiply_int; + math->md_divide_int = mp_number_divide_int; + math->md_to_boolean = mp_number_to_boolean; + math->md_to_scaled = mp_number_to_scaled; + math->md_to_double = mp_number_to_double; + math->md_to_int = mp_number_to_int; + math->md_odd = mp_number_odd; + math->md_equal = mp_number_equal; + math->md_less = mp_number_less; + math->md_greater = mp_number_greater; + math->md_nonequalabs = mp_number_nonequalabs; + math->md_round_unscaled = mp_round_unscaled; + math->md_floor_scaled = mp_number_floor; + math->md_fraction_to_round_scaled = mp_posit_fraction_to_round_scaled; + math->md_make_scaled = mp_posit_number_make_scaled; + math->md_make_fraction = mp_posit_number_make_fraction; + math->md_take_fraction = mp_posit_number_take_fraction; + math->md_take_scaled = mp_posit_number_take_scaled; + math->md_velocity = mp_posit_velocity; + math->md_n_arg = mp_posit_n_arg; + math->md_m_log = mp_posit_m_log; + math->md_m_exp = mp_posit_m_exp; + math->md_m_unif_rand = mp_posit_m_unif_rand; + math->md_m_norm_rand = mp_posit_m_norm_rand; + math->md_pyth_add = mp_posit_pyth_add; + math->md_pyth_sub = mp_posit_pyth_sub; + math->md_power_of = mp_posit_power_of; + math->md_fraction_to_scaled = mp_number_fraction_to_scaled; + math->md_scaled_to_fraction = mp_number_scaled_to_fraction; + math->md_scaled_to_angle = mp_number_scaled_to_angle; + math->md_angle_to_scaled = mp_number_angle_to_scaled; + math->md_init_randoms = mp_init_randoms; + math->md_sin_cos = mp_posit_sin_cos; + math->md_slow_add = mp_posit_slow_add; + math->md_sqrt = mp_posit_square_rt; + math->md_print = mp_posit_print_number; + math->md_tostring = mp_posit_number_tostring; + math->md_modulo = mp_number_modulo; + math->md_ab_vs_cd = mp_posit_ab_vs_cd; + math->md_crossing_point = mp_posit_crossing_point; + math->md_scan_numeric = mp_posit_scan_numeric_token; + math->md_scan_fractional = mp_posit_scan_fractional_token; + math->md_free_math = mp_free_posit_math; + math->md_set_precision = mp_posit_set_precision; + return math; +} + +void mp_posit_set_precision (MP mp) +{ + (void) mp; +} + +void mp_free_posit_math (MP mp) +{ + /* Is this list up to date? Also check elewhere. */ + mp_free_number(mp, &(mp->math->md_three_sixty_deg_t)); + mp_free_number(mp, &(mp->math->md_one_eighty_deg_t)); + mp_free_number(mp, &(mp->math->md_negative_one_eighty_deg_t)); + mp_free_number(mp, &(mp->math->md_fraction_one_t)); + mp_free_number(mp, &(mp->math->md_zero_t)); + mp_free_number(mp, &(mp->math->md_half_unit_t)); + mp_free_number(mp, &(mp->math->md_three_quarter_unit_t)); + mp_free_number(mp, &(mp->math->md_unity_t)); + mp_free_number(mp, &(mp->math->md_two_t)); + mp_free_number(mp, &(mp->math->md_three_t)); + mp_free_number(mp, &(mp->math->md_one_third_inf_t)); + mp_free_number(mp, &(mp->math->md_inf_t)); + mp_free_number(mp, &(mp->math->md_negative_inf_t)); + mp_free_number(mp, &(mp->math->md_warning_limit_t)); + mp_free_number(mp, &(mp->math->md_one_k)); + mp_free_number(mp, &(mp->math->md_sqrt_8_e_k)); + mp_free_number(mp, &(mp->math->md_twelve_ln_2_k)); + mp_free_number(mp, &(mp->math->md_coef_bound_k)); + mp_free_number(mp, &(mp->math->md_coef_bound_minus_1)); + mp_free_number(mp, &(mp->math->md_fraction_threshold_t)); + mp_free_number(mp, &(mp->math->md_half_fraction_threshold_t)); + mp_free_number(mp, &(mp->math->md_scaled_threshold_t)); + mp_free_number(mp, &(mp->math->md_half_scaled_threshold_t)); + mp_free_number(mp, &(mp->math->md_near_zero_angle_t)); + mp_free_number(mp, &(mp->math->md_p_over_v_threshold_t)); + mp_free_number(mp, &(mp->math->md_equation_threshold_t)); + mp_memory_free(mp->math); +} + +@ See mpmathdouble for documentation. @c + +void mp_allocate_number (MP mp, mp_number *n, mp_number_type t) +{ + (void) mp; + n->data.pval = mp_posit_data.zero; + n->type = t; +} + +void mp_allocate_clone (MP mp, mp_number *n, mp_number_type t, mp_number *v) +{ + (void) mp; + n->type = t; + n->data.pval = v->data.pval; +} + +void mp_allocate_abs (MP mp, mp_number *n, mp_number_type t, mp_number *v) +{ + (void) mp; + n->type = t; + n->data.pval = posit_fabs(v->data.pval); +} + +void mp_allocate_double (MP mp, mp_number *n, double v) +{ + (void) mp; + n->type = mp_scaled_type; + n->data.pval = double_to_posit(v); +} + +void mp_free_number (MP mp, mp_number *n) +{ + (void) mp; + n->type = mp_nan_type; +} + +void mp_set_posit_from_int(mp_number *A, int B) +{ + A->data.pval = integer_to_posit(B); +} + +void mp_set_posit_from_boolean(mp_number *A, int B) +{ + A->data.pval = integer_to_posit(B); +} + +void mp_set_posit_from_scaled(mp_number *A, int B) +{ + A->data.pval = posit_div(integer_to_posit(B), mp_posit_data.d65536); +} + +void mp_set_posit_from_double(mp_number *A, double B) +{ + A->data.pval = double_to_posit(B); +} + +void mp_set_posit_from_addition(mp_number *A, mp_number *B, mp_number *C) +{ + A->data.pval = posit_add(B->data.pval, C->data.pval); +} + +void mp_set_posit_half_from_addition(mp_number *A, mp_number *B, mp_number *C) +{ + A->data.pval = posit_div(posit_add(B->data.pval,C->data.pval), mp_posit_data.two); +} + +void mp_set_posit_from_subtraction(mp_number *A, mp_number *B, mp_number *C) +{ + A->data.pval = posit_sub(B->data.pval, C->data.pval); +} + +void mp_set_posit_half_from_subtraction(mp_number *A, mp_number *B, mp_number *C) +{ + A->data.pval = posit_div(posit_sub(B->data.pval, C->data.pval), mp_posit_data.two); +} + +void mp_set_posit_from_div(mp_number *A, mp_number *B, mp_number *C) +{ + A->data.pval = posit_div(B->data.pval, C->data.pval); +} + +void mp_set_posit_from_mul(mp_number *A, mp_number *B, mp_number *C) +{ + A->data.pval = posit_mul(B->data.pval, C->data.pval); +} + +void mp_set_posit_from_int_div(mp_number *A, mp_number *B, int C) +{ + A->data.pval = posit_div(B->data.pval, integer_to_posit(C)); +} + +void mp_set_posit_from_int_mul(mp_number *A, mp_number *B, int C) +{ + A->data.pval = posit_mul(A->data.pval, integer_to_posit(C)); +} + +void mp_set_posit_from_of_the_way (MP mp, mp_number *A, mp_number *t, mp_number *B, mp_number *C) +{ + (void) mp; + A->data.pval = posit_sub(B->data.pval, mp_posit_take_fraction(posit_sub(B->data.pval, C->data.pval), t->data.pval)); +} + +void mp_number_negate(mp_number *A) +{ + A->data.pval = posit_neg(A->data.pval); +} + +void mp_number_add(mp_number *A, mp_number *B) +{ + A->data.pval = posit_add(A->data.pval, B->data.pval); +} + +void mp_number_subtract(mp_number *A, mp_number *B) +{ + A->data.pval = posit_sub(A->data.pval, B->data.pval); +} + +void mp_number_half(mp_number *A) +{ + A->data.pval = posit_div(A->data.pval, mp_posit_data.two); +} + +void mp_number_double(mp_number *A) +{ + A->data.pval = posit_mul(A->data.pval, mp_posit_data.two); +} + +void mp_number_add_scaled(mp_number *A, int B) +{ + /* also for negative B */ + A->data.pval = posit_add(A->data.pval, posit_div(integer_to_posit(B), mp_posit_data.d65536)); +} + +void mp_number_multiply_int(mp_number *A, int B) +{ + A->data.pval = posit_mul(A->data.pval, integer_to_posit(B)); +} + +void mp_number_divide_int(mp_number *A, int B) +{ + A->data.pval = posit_div(A->data.pval, integer_to_posit(B)); +} + +void mp_posit_abs(mp_number *A) +{ + A->data.pval = posit_fabs(A->data.pval); +} + +void mp_number_clone(mp_number *A, mp_number *B) +{ + A->data.pval = B->data.pval; +} + +void mp_number_negated_clone(mp_number *A, mp_number *B) +{ + A->data.pval = posit_neg(B->data.pval); +} + +void mp_number_abs_clone(mp_number *A, mp_number *B) +{ + A->data.pval = posit_fabs(B->data.pval); +} + +void mp_number_swap(mp_number *A, mp_number *B) +{ + posit_t swap_tmp = A->data.pval; + A->data.pval = B->data.pval; + B->data.pval = swap_tmp; +} + +void mp_number_fraction_to_scaled(mp_number *A) +{ + A->type = mp_scaled_type; + A->data.pval = posit_div(A->data.pval, mp_posit_data.fraction_multiplier); +} + +void mp_number_angle_to_scaled(mp_number *A) +{ + A->type = mp_scaled_type; + A->data.pval = posit_div(A->data.pval, mp_posit_data.angle_multiplier); +} + +void mp_number_scaled_to_fraction(mp_number *A) +{ + A->type = mp_fraction_type; + A->data.pval = posit_mul(A->data.pval, mp_posit_data.fraction_multiplier); +} + +void mp_number_scaled_to_angle(mp_number *A) +{ + A->type = mp_angle_type; + A->data.pval = posit_mul(A->data.pval, mp_posit_data.angle_multiplier); +} + +int mp_number_to_scaled(mp_number *A) +{ + return posit_to_integer(posit_mul(A->data.pval, mp_posit_data.d65536)); +} + +int mp_number_to_int(mp_number *A) +{ + return posit_to_integer(A->data.pval); +} + +int mp_number_to_boolean(mp_number *A) +{ + return posit_eq_zero(A->data.pval) ? 0 : 1; +} + +double mp_number_to_double(mp_number *A) +{ + return posit_to_double(A->data.pval); +} + +int mp_number_odd(mp_number *A) +{ + return odd(posit_to_integer(A->data.pval)); +} + +int mp_number_equal(mp_number *A, mp_number *B) +{ + return posit_eq(A->data.pval, B->data.pval); +} + +int mp_number_greater(mp_number *A, mp_number *B) +{ + return posit_gt(A->data.pval, B->data.pval); +} + +int mp_number_less(mp_number *A, mp_number *B) +{ + return posit_lt(A->data.pval, B->data.pval); +} + +int mp_number_nonequalabs(mp_number *A, mp_number *B) +{ + return ! posit_eq(posit_fabs(A->data.pval), posit_fabs(B->data.pval)); +} + +char *mp_posit_number_tostring (MP mp, mp_number *n) +{ + static char set[64]; + int l = 0; + char *ret = mp_memory_allocate(64); + (void) mp; + snprintf(set, 64, "%.20g", posit_to_double(n->data.pval)); + while (set[l] == ' ') { + l++; + } + strcpy(ret, set+l); + return ret; +} + +void mp_posit_print_number (MP mp, mp_number *n) +{ + char *str = mp_posit_number_tostring(mp, n); + mp_print_e_str(mp, str); + mp_memory_free(str); +} + +/* Todo: it is hard to overflow posits. Also, we can check zero fast. */ + +void mp_posit_slow_add (MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_orig) +{ + if (posit_gt(x_orig->data.pval, mp_posit_data.zero)) { + if (posit_le(y_orig->data.pval, posit_sub(mp_posit_data.EL_GORDO, x_orig->data.pval))) { + ret->data.pval = posit_add(x_orig->data.pval, y_orig->data.pval); + } else { + mp->arith_error = 1; + ret->data.pval = mp_posit_data.EL_GORDO; + } + } else if (posit_le(posit_neg(y_orig->data.pval), posit_add(mp_posit_data.EL_GORDO, x_orig->data.pval))) { + ret->data.pval = posit_add(x_orig->data.pval, y_orig->data.pval); + } else { + mp->arith_error = 1; + ret->data.pval = mp_posit_data.negative_EL_GORDO; + } +} + +void mp_posit_number_make_fraction (MP mp, mp_number *ret, mp_number *p, mp_number *q) { + (void) mp; + ret->data.pval = mp_posit_make_fraction(p->data.pval, q->data.pval); +} + +void mp_posit_number_take_fraction (MP mp, mp_number *ret, mp_number *p, mp_number *q) { + (void) mp; + ret->data.pval = mp_posit_take_fraction(p->data.pval, q->data.pval); +} + +void mp_posit_number_take_scaled (MP mp, mp_number *ret, mp_number *p_orig, mp_number *q_orig) +{ + (void) mp; + ret->data.pval = posit_mul(p_orig->data.pval, q_orig->data.pval); +} + +void mp_posit_number_make_scaled (MP mp, mp_number *ret, mp_number *p_orig, mp_number *q_orig) +{ + (void) mp; + ret->data.pval = posit_div(p_orig->data.pval, q_orig->data.pval); +} + +void mp_wrapup_numeric_token (MP mp, unsigned char *start, unsigned char *stop) +{ + double result; + char *end = (char *) stop; + errno = 0; + result = strtod((char *) start, &end); + if (errno == 0) { + set_cur_mod(double_to_posit(result)); + if (result >= mp_warning_limit) { + if (posit_gt(internal_value(mp_warning_check_internal).data.pval, mp_posit_data.zero) && (mp->scanner_status != mp_tex_flushing_state)) { + char msg[256]; + mp_snprintf(msg, 256, "Number is too large (%g)", result); + @.Number is too large@> + mp_error( + mp, + msg, + "Continue and I'll try to cope with that big value; but it might be dangerous." + "(Set warningcheck := 0 to suppress this message.)" + ); + } + } + } else if (mp->scanner_status != mp_tex_flushing_state) { + mp_error( + mp, + "Enormous number has been reduced.", + "I could not handle this number specification probably because it is out of" + "range." + ); + @.Enormous number...@> + set_cur_mod(mp_posit_data.EL_GORDO); + } + set_cur_cmd(mp_numeric_command); +} + +static void mp_posit_aux_find_exponent (MP mp) +{ + if (mp->buffer[mp->cur_input.loc_field] == 'e' || mp->buffer[mp->cur_input.loc_field] == 'E') { + mp->cur_input.loc_field++; + if (!(mp->buffer[mp->cur_input.loc_field] == '+' + || mp->buffer[mp->cur_input.loc_field] == '-' + || mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class)) { + mp->cur_input.loc_field--; + return; + } + if (mp->buffer[mp->cur_input.loc_field] == '+' + || mp->buffer[mp->cur_input.loc_field] == '-') { + mp->cur_input.loc_field++; + } + while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class) { + mp->cur_input.loc_field++; + } + } +} + +void mp_posit_scan_fractional_token (MP mp, int n) /* n is scaled */ +{ + unsigned char *start = &mp->buffer[mp->cur_input.loc_field -1]; + unsigned char *stop; + (void) n; + while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class) { + mp->cur_input.loc_field++; + } + mp_posit_aux_find_exponent(mp); + stop = &mp->buffer[mp->cur_input.loc_field-1]; + mp_wrapup_numeric_token(mp, start, stop); +} + +void mp_posit_scan_numeric_token (MP mp, int n) /* n is scaled */ +{ + unsigned char *start = &mp->buffer[mp->cur_input.loc_field -1]; + unsigned char *stop; + (void) n; + while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class) { + mp->cur_input.loc_field++; + } + if (mp->buffer[mp->cur_input.loc_field] == '.' && mp->buffer[mp->cur_input.loc_field+1] != '.') { + mp->cur_input.loc_field++; + while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class) { + mp->cur_input.loc_field++; + } + } + mp_posit_aux_find_exponent(mp); + stop = &mp->buffer[mp->cur_input.loc_field-1]; + mp_wrapup_numeric_token(mp, start, stop); +} + +void mp_posit_velocity (MP mp, mp_number *ret, mp_number *st, mp_number *ct, mp_number *sf, mp_number *cf, mp_number *t) +{ + posit_t acc, num, denom; /* registers for intermediate calculations */ + (void) mp; + acc = mp_posit_take_fraction( + mp_posit_take_fraction( + posit_sub(st->data.pval, posit_div(sf->data.pval, mp_posit_data.sixteen)), + posit_sub(sf->data.pval, posit_div(st->data.pval, mp_posit_data.sixteen)) + ), + posit_sub(ct->data.pval,cf->data.pval) + ); + num = posit_add( + mp_posit_data.fraction_two, + mp_posit_take_fraction( + acc, + mp_posit_data.sqrt_two_mul_fraction_one + ) + ); + denom = posit_add( + mp_posit_data.fraction_three, + posit_add( + mp_posit_take_fraction( + ct->data.pval, + mp_posit_data.sqrt_five_minus_one_mul_fraction_one_and_half + ), + mp_posit_take_fraction( + cf->data.pval, + mp_posit_data.three_minus_sqrt_five_mul_fraction_one_and_half + ) + ) + ); + if (posit_ne(t->data.pval, mp_posit_data.unity)) { + num = mp_posit_make_scaled(num, t->data.pval); + } + if (posit_ge(posit_div(num, mp_posit_data.four), denom)) { + ret->data.pval = mp_posit_data.fraction_four; + } else { + ret->data.pval = mp_posit_make_fraction(num, denom); + } +} + +int mp_posit_ab_vs_cd (mp_number *a_orig, mp_number *b_orig, mp_number *c_orig, mp_number *d_orig) +{ + posit_t ab = posit_mul(a_orig->data.pval, b_orig->data.pval); + posit_t cd = posit_mul(c_orig->data.pval, d_orig->data.pval); + if (posit_eq(ab,cd)) { + return 0; + } else if (posit_lt(ab,cd)) { + return -1; + } else { + return 1; + } +} + +static void mp_posit_crossing_point (MP mp, mp_number *ret, mp_number *aa, mp_number *bb, mp_number *cc) +{ + posit_t d; + posit_t xx, x0, x1, x2; + posit_t a = aa->data.pval; + posit_t b = bb->data.pval; + posit_t c = cc->data.pval; + (void) mp; + if (posit_lt(a, mp_posit_data.zero)) { + ret->data.pval = mp_posit_data.zero_crossing; + return; + } + if (posit_ge(c, mp_posit_data.zero)) { + if (posit_ge(b, mp_posit_data.zero)) { + if (posit_gt(c, mp_posit_data.zero)) { + ret->data.pval = mp_posit_data.no_crossing; + } else if (posit_eq_zero(a) && posit_eq_zero(b)) { + ret->data.pval = mp_posit_data.no_crossing; + } else { + ret->data.pval = mp_posit_data.one_crossing; + } + return; + } + if (posit_eq_zero(a)) { + ret->data.pval = mp_posit_data.zero_crossing; + return; + } + } else if (posit_eq_zero(a) && posit_le(b, mp_posit_data.zero)) { + ret->data.pval = mp_posit_data.zero_crossing; + return; + } + /* Use bisection to find the crossing point... */ + d = mp_posit_data.epsilon; + x0 = a; + x1 = posit_sub(a, b); + x2 = posit_sub(b, c); + do { + /* not sure why the error correction has to be >= 1E-12 */ + posit_t x = posit_add(posit_div(posit_add(x1, x2), mp_posit_data.two), mp_posit_data.error_correction); + if (posit_gt(posit_sub(x1, x0), x0)) { + x2 = x; + x0 = posit_add(x0, x0); + d = posit_add(d, d); + } else { + xx = posit_sub(posit_add(x1, x), x0); + if (posit_gt(xx, x0)) { + x2 = x; + x0 = posit_add(x0, x0); + d = posit_add(d, d); + } else { + x0 = posit_sub(x0, xx); + if (posit_le(x, x0) && posit_le(posit_add(x, x2), x0)) { + ret->data.pval = mp_posit_data.no_crossing; + return; + } + x1 = x; + d = posit_add(posit_add(d, d), mp_posit_data.epsilon); + } + } + } while (posit_lt(d, mp_posit_data.fraction_one)); + ret->data.pval = posit_sub(d, mp_posit_data.fraction_one); +} + +@ See mpmathdouble for documentation. @c + +int mp_round_unscaled(mp_number *x_orig) +{ + return posit_i_round(x_orig->data.pval); +} + +void mp_number_floor(mp_number *i) +{ + i->data.pval = posit_floor(i->data.pval); +} + +void mp_posit_fraction_to_round_scaled(mp_number *x_orig) +{ + x_orig->type = mp_scaled_type; + x_orig->data.pval = posit_div(x_orig->data.pval, mp_posit_data.fraction_multiplier); +} + +void mp_posit_square_rt (MP mp, mp_number *ret, mp_number *x_orig) /* return, x: scaled */ +{ + if (posit_gt(x_orig->data.pval, mp_posit_data.zero)) { + ret->data.pval = posit_sqrt(x_orig->data.pval); + } else { + if (posit_lt(x_orig->data.pval, mp_posit_data.zero)) { + char msg[256]; + char *xstr = mp_posit_number_tostring(mp, x_orig); + mp_snprintf(msg, 256, "Square root of %s has been replaced by 0", xstr); + mp_memory_free(xstr); + @.Square root...replaced by 0@> + mp_error( + mp, + msg, + "Since I don't take square roots of negative numbers, I'm zeroing this one.\n" + "Proceed, with fingers crossed." + ); + } + ret->data.pval = mp_posit_data.zero; + } +} + +void mp_posit_pyth_add (MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig) +{ + ret->data.pval = posit_sqrt( + posit_add( + posit_mul( + a_orig->data.pval, + a_orig->data.pval + ), + posit_mul( + b_orig->data.pval, + b_orig->data.pval + ) + ) + ); +} + +void mp_posit_pyth_sub (MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig) +{ + /* can be made nicer */ + if (posit_gt(a_orig->data.pval,b_orig->data.pval)) { + a_orig->data.pval = posit_sqrt( + posit_sub( + posit_mul( + a_orig->data.pval, + a_orig->data.pval + ), + posit_mul( + b_orig->data.pval, + b_orig->data.pval + ) + ) + ); + } else { + if (posit_lt(a_orig->data.pval,b_orig->data.pval)) { + char msg[256]; + char *astr = mp_posit_number_tostring(mp, a_orig); + char *bstr = mp_posit_number_tostring(mp, b_orig); + mp_snprintf(msg, 256, "Pythagorean subtraction %s+-+%s has been replaced by 0", astr, bstr); + mp_memory_free(astr); + mp_memory_free(bstr); + @.Pythagorean...@> + mp_error( + mp, + msg, + "Since I don't take square roots of negative numbers, Im zeroing this one.\n" + "Proceed, with fingers crossed." + ); + } + a_orig->data.pval = mp_posit_data.zero; + } + ret->data.pval = a_orig->data.pval; +} + +void mp_posit_power_of (MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig) +{ + errno = 0; + ret->data.pval = posit_pow(a_orig->data.pval, b_orig->data.pval); + if (errno) { + mp->arith_error = 1; + ret->data.pval = mp_posit_data.EL_GORDO; + } +} + +void mp_posit_m_log (MP mp, mp_number *ret, mp_number *x_orig) +{ + /* TODO: int mult */ + if (posit_gt(x_orig->data.pval,mp_posit_data.zero)) { + ret->data.pval = posit_mul(posit_log(x_orig->data.pval),mp_posit_data.d256); + } else { + char msg[256]; + char *xstr = mp_posit_number_tostring(mp, x_orig); + mp_snprintf(msg, 256, "Logarithm of %s has been replaced by 0", xstr); + mp_memory_free(xstr); + mp_error( + mp, + msg, + "Since I don't take logs of non-positive numbers, I'm zeroing this one.\n" + "Proceed, with fingers crossed." + ); + ret->data.pval = mp_posit_data.zero; + } +} + +void mp_posit_m_exp (MP mp, mp_number *ret, mp_number *x_orig) +{ + errno = 0; + ret->data.pval = posit_exp(posit_div(x_orig->data.pval,mp_posit_data.d256)); + if (errno) { + if (posit_gt(x_orig->data.pval,mp_posit_data.zero)) { + mp->arith_error = 1; + ret->data.pval = mp_posit_data.EL_GORDO; + } else { + ret->data.pval = mp_posit_data.zero; + } + } +} + +void mp_posit_n_arg (MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_orig) +{ + if (posit_eq_zero(x_orig->data.pval) && posit_eq_zero(y_orig->data.pval)) { + mp_error( + mp, + "angle(0,0) is taken as zero", + "The 'angle' between two identical points is undefined. I'm zeroing this one.\n" + "Proceed, with fingers crossed." + ); + ret->data.pval = mp_posit_data.zero; + } else { + ret->type = mp_angle_type; + /* TODO */ + ret->data.pval = posit_mul( + posit_atan2( + y_orig->data.pval, + x_orig->data.pval + ), + mp_posit_data.d180_divided_by_pi_mul_angle + ); + } +} + +void mp_posit_sin_cos (MP mp, mp_number *z_orig, mp_number *n_cos, mp_number *n_sin) +{ + posit_t rad = posit_div(z_orig->data.pval, mp_posit_data.angle_multiplier); + (void) mp; + if (posit_eq(rad, mp_posit_data.dp90) || posit_eq(rad, mp_posit_data.dm270)) { + n_cos->data.pval = mp_posit_data.zero; + n_sin->data.pval = mp_posit_data.fraction_multiplier; + } else if (posit_eq(rad, mp_posit_data.dm90) || posit_eq(rad, mp_posit_data.dp270)) { + n_cos->data.pval = mp_posit_data.zero; + n_sin->data.pval = mp_posit_data.negative_fraction_multiplier; + } else if (posit_eq(rad, mp_posit_data.dp180) || posit_eq(rad, mp_posit_data.dm180)) { + n_cos->data.pval = mp_posit_data.negative_fraction_multiplier; + n_sin->data.pval = mp_posit_data.zero; + } else { + rad = posit_mul(rad,mp_posit_data.pi_divided_by_180); + n_cos->data.pval = posit_mul(posit_cos(rad),mp_posit_data.fraction_multiplier); + n_sin->data.pval = posit_mul(posit_sin(rad),mp_posit_data.fraction_multiplier); + } +} + +@ See mpmathdouble for documentation. @c + +# define KK 100 /* the long lag */ +# define LL 37 /* the short lag */ +# define MM (1L<<30) /* the modulus */ +# define mod_diff(x,y) (((x)-(y))&(MM-1)) /* subtraction mod MM */ +# define TT 70 /* guaranteed separation between streams */ +# define is_odd(x) ((x)&1) /* units bit of x */ +# define QUALITY 1009 /* recommended quality level for high-res use */ + +/* destination, array length (must be at least KK) */ + +typedef struct mp_posit_random_info { + long x[KK]; + long buf[QUALITY]; + long dummy; + long started; + long *ptr; +} mp_posit_random_info; + +static mp_posit_random_info mp_posit_random_data = { + .dummy = -1, + .started = -1, + .ptr = &mp_posit_random_data.dummy +}; + +/* the following routines are from exercise 3.6--15 */ +/* after calling |mp_aux_ran_start|, get new randoms by, e.g., |x=mp_aux_ran_arr_next()| */ + +static void mp_posit_aux_ran_array(long aa[], int n) +{ + int i, j; + for (j = 0; j < KK; j++) { + aa[j] = mp_posit_random_data.x[j]; + } + for (; j < n; j++) { + aa[j] = mod_diff(aa[j - KK], aa[j - LL]); + } + for (i = 0; i < LL; i++, j++) { + mp_posit_random_data.x[i] = mod_diff(aa[j - KK], aa[j - LL]); + } + for (; i < KK; i++, j++) { + mp_posit_random_data.x[i] = mod_diff(aa[j - KK], mp_posit_random_data.x[i - LL]); + } +} + +/* Do this before using |mp_aux_ran_array|, long seed selector for different streams. */ + +static void mp_posit_aux_ran_start(long seed) +{ + int t, j; + long x[KK + KK - 1]; /* the preparation buffer */ + long ss = (seed+2) & (MM - 2); + for (j = 0; j < KK; j++) { + /* bootstrap the buffer */ + x[j] = ss; + /* cyclic shift 29 bits */ + ss <<= 1; + if (ss >= MM) { + ss -= MM - 2; + } + } + /* make x[1] (and only x[1]) odd */ + x[1]++; + for (ss = seed & (MM - 1), t = TT - 1; t;) { + for (j = KK - 1; j > 0; j--) { + /* "square" */ + x[j + j] = x[j]; + x[j + j - 1] = 0; + } + for (j = KK + KK - 2; j >= KK; j--) { + x[j - (KK -LL)] = mod_diff(x[j - (KK - LL)], x[j]); + x[j - KK] = mod_diff(x[j - KK], x[j]); + } + if (is_odd(ss)) { + /* "multiply by z" */ + for (j = KK; j>0; j--) { + x[j] = x[j-1]; + } + x[0] = x[KK]; + /* shift the buffer cyclically */ + x[LL] = mod_diff(x[LL], x[KK]); + } + if (ss) { + ss >>= 1; + } else { + t--; + } + } + for (j = 0; j < LL; j++) { + mp_posit_random_data.x[j + KK - LL] = x[j]; + } + for (;j < KK; j++) { + mp_posit_random_data.x[j - LL] = x[j]; + } + for (j = 0; j < 10; j++) { + /* warm things up */ + mp_posit_aux_ran_array(x, KK + KK - 1); + } + mp_posit_random_data.ptr = &mp_posit_random_data.started; +} + +static long mp_posit_aux_ran_arr_cycle(void) +{ + if (mp_posit_random_data.ptr == &mp_posit_random_data.dummy) { + /* the user forgot to initialize */ + mp_posit_aux_ran_start(314159L); + } + mp_posit_aux_ran_array(mp_posit_random_data.buf, QUALITY); + mp_posit_random_data.buf[KK] = -1; + mp_posit_random_data.ptr = mp_posit_random_data.buf + 1; + return mp_posit_random_data.buf[0]; +} + +void mp_init_randoms (MP mp, int seed) +{ + int k = 1; + int j = abs(seed); + int f = (int) mp_fraction_multiplier; /* avoid warnings */ + while (j >= f) { + j = j/2; + } + for (int i = 0; i <= 54; i++) { + int jj = k; + k = j - k; + j = jj; + if (k < 0) { + k += f; + } + mp->randoms[(i * 21) % 55].data.pval = integer_to_posit(j); + } + mp_new_randoms(mp); + mp_new_randoms(mp); + mp_new_randoms(mp); + /* warm up the array */ + mp_posit_aux_ran_start((unsigned long) seed); +} + +void mp_number_modulo(mp_number *a, mp_number *b) +{ + a->data.pval = posit_mul(posit_modf(posit_div(a->data.pval, b->data.pval)), b->data.pval); +} + +static void mp_next_unif_random (MP mp, mp_number *ret) +{ + unsigned long int op = (unsigned) (*mp_posit_random_data.ptr >=0 ? *mp_posit_random_data.ptr++: mp_posit_aux_ran_arr_cycle()); + double a = op / (MM * 1.0); + (void) mp; + ret->data.pval = double_to_posit(a); +} + +static void mp_next_random (MP mp, mp_number *ret) +{ + if ( mp->j_random==0) { + mp_new_randoms(mp); + } else { + mp->j_random = mp->j_random-1; + } + mp_number_clone(ret, &(mp->randoms[mp->j_random])); +} + +static void mp_posit_m_unif_rand (MP mp, mp_number *ret, mp_number *x_orig) +{ + mp_number x, abs_x, u, y; + mp_allocate_number(mp, &y, mp_fraction_type); + mp_allocate_clone(mp, &x, mp_scaled_type, x_orig); + mp_allocate_abs(mp, &abs_x, mp_scaled_type, &x); + mp_allocate_number(mp, &u, mp_scaled_type); + mp_next_unif_random(mp, &u); + y.data.pval = posit_mul(abs_x.data.pval, u.data.pval); + mp_free_number(mp, &u); + if (mp_number_equal(&y, &abs_x)) { + mp_number_clone(ret, &((math_data *)mp->math)->md_zero_t); + } else if (mp_number_greater(&x, &((math_data *)mp->math)->md_zero_t)) { + mp_number_clone(ret, &y); + } else { + mp_number_negated_clone(ret, &y); + } + mp_free_number(mp, &abs_x); + mp_free_number(mp, &x); + mp_free_number(mp, &y); +} + +static void mp_posit_m_norm_rand (MP mp, mp_number *ret) +{ + mp_number abs_x, u, r, la, xa; + mp_allocate_number(mp, &la, mp_scaled_type); + mp_allocate_number(mp, &xa, mp_scaled_type); + mp_allocate_number(mp, &abs_x, mp_scaled_type); + mp_allocate_number(mp, &u, mp_scaled_type); + mp_allocate_number(mp, &r, mp_scaled_type); + do { + do { + mp_number v; + mp_allocate_number(mp, &v, mp_scaled_type); + mp_next_random(mp, &v); + mp_number_subtract(&v, &((math_data *)mp->math)->md_fraction_half_t); + mp_posit_number_take_fraction(mp, &xa, &((math_data *)mp->math)->md_sqrt_8_e_k, &v); + mp_free_number(mp, &v); + mp_next_random(mp, &u); + mp_number_clone(&abs_x, &xa); + mp_posit_abs(&abs_x); + } while (! mp_number_less(&abs_x, &u)); + mp_posit_number_make_fraction(mp, &r, &xa, &u); + mp_number_clone(&xa, &r); + mp_posit_m_log(mp, &la, &u); + mp_set_posit_from_subtraction(&la, &((math_data *)mp->math)->md_twelve_ln_2_k, &la); + } while (mp_posit_ab_vs_cd(&((math_data *)mp->math)->md_one_k, &la, &xa, &xa) < 0); + mp_number_clone(ret, &xa); + mp_free_number(mp, &r); + mp_free_number(mp, &abs_x); + mp_free_number(mp, &la); + mp_free_number(mp, &xa); + mp_free_number(mp, &u); +} + +@ @<Reduce to the case that |a...@>= +if (posit_lt(a, mp_posit_data.zero)) { + a = posit_neg(a); + b = posit_neg(b); +} +if (posit_lt(c, mp_posit_data.zero)) { + c = posit_neg(c); + d = posit_neg(d); +} +if ((posit_le(d, mp_posit_data.zero)) { + if ((posit_ge(b, mp_posit_data.zero)) { + if ((posit_eq_zero(a) || posit_eq_zero(b) && (posit_eq_zero(c) || posit_eq_zero(d))) { + ret->data.pval = mp_posit_data.zero; + } else { + ret->data.pval = mp_posit_data.one; + } + goto RETURN; + } if (posit_eq_zero(d)) { + ret->data.pval = posit_eq_zero(a) ? mp_posit_data.zero : mp_posit_data.minusone; + goto RETURN; + } else + q = a; + a = c; + c = q; + q = -b; + b = -d; + d = q; + } +} else if (posit_le(b, mp_posit_data.zero) { + if (posit_lt(b, mp_posit_data.zero) && posit_gt(a, mp_posit_data.zero)) { + ret->data.pval = mp_posit_data.minusone; + return; + } else + ret->data.pval = posit_eq_zero(c) ? mp_posit_data.zero : mp_posit_data.minusone; + goto RETURN; + } +} |