?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
beta-0.89.2
[luatex.git] / source / texk / web2c / mplibdir / mpmathbinary.w
blob4d7301579a1c77cd54eb4f67bd740ebb12094eee
1 % $Id$
2 %
3 % This file is part of MetaPost;
4 % the MetaPost program is in the public domain.
5 % See the <Show version...> code in mpost.w for more info.
6
7 \def\title{Math support functions for MPFR based math}
8 \pdfoutput=1
9
10 @ Introduction.
11
12 @c
13 #include <w2c/config.h>
14 #include <stdio.h>
15 #include <stdlib.h>
16 #include <string.h>
17 #include <math.h>
18 #include "mpmathbinary.h" /* internal header */
19 #define ROUND(a) floor((a)+0.5)
20 @h
21
22 @ @c
23 @<Declarations@>;
24
25 @ @(mpmathbinary.h@>=
26 #ifndef MPMATHBINARY_H
27 #define MPMATHBINARY_H 1
28 #include "mplib.h"
29 #include "mpmp.h" /* internal header */
30 #include <gmp.h>
31 #include <mpfr.h>
32 @<Internal library declarations@>;
33 #endif
34
35 @* Math initialization.
36
37 First, here are some very important constants.
38
39 @d ROUNDING MPFR_RNDN
40 @d E_STRING "2.7182818284590452353602874713526624977572470936999595749669676277240766303535"
41 @d PI_STRING "3.1415926535897932384626433832795028841971693993751058209749445923078164062862"
42 @d fraction_multiplier 4096
43 @d angle_multiplier 16
44
45 @ Here are the functions that are static as they are not used elsewhere
46
47 @<Declarations@>=
48 #define DEBUG 0
49 static void mp_binary_scan_fractional_token (MP mp, int n);
50 static void mp_binary_scan_numeric_token (MP mp, int n);
51 static void mp_binary_ab_vs_cd (MP mp, mp_number *ret, mp_number a, mp_number b, mp_number c, mp_number d);
52 static void mp_ab_vs_cd (MP mp, mp_number *ret, mp_number a, mp_number b, mp_number c, mp_number d);
53 static void mp_binary_crossing_point (MP mp, mp_number *ret, mp_number a, mp_number b, mp_number c);
54 static void mp_binary_number_modulo (mp_number *a, mp_number b);
55 static void mp_binary_print_number (MP mp, mp_number n);
56 static char * mp_binary_number_tostring (MP mp, mp_number n);
57 static void mp_binary_slow_add (MP mp, mp_number *ret, mp_number x_orig, mp_number y_orig);
58 static void mp_binary_square_rt (MP mp, mp_number *ret, mp_number x_orig);
59 static void mp_binary_sin_cos (MP mp, mp_number z_orig, mp_number *n_cos, mp_number *n_sin);
60 static void mp_init_randoms (MP mp, int seed);
61 static void mp_number_angle_to_scaled (mp_number *A);
62 static void mp_number_fraction_to_scaled (mp_number *A);
63 static void mp_number_scaled_to_fraction (mp_number *A);
64 static void mp_number_scaled_to_angle (mp_number *A);
65 static void mp_binary_m_unif_rand (MP mp, mp_number *ret, mp_number x_orig);
66 static void mp_binary_m_norm_rand (MP mp, mp_number *ret);
67 static void mp_binary_m_exp (MP mp, mp_number *ret, mp_number x_orig);
68 static void mp_binary_m_log (MP mp, mp_number *ret, mp_number x_orig);
69 static void mp_binary_pyth_sub (MP mp, mp_number *r, mp_number a, mp_number b);
70 static void mp_binary_pyth_add (MP mp, mp_number *r, mp_number a, mp_number b);
71 static void mp_binary_n_arg (MP mp, mp_number *ret, mp_number x, mp_number y);
72 static void mp_binary_velocity (MP mp, mp_number *ret, mp_number st, mp_number ct, mp_number sf, mp_number cf, mp_number t);
73 static void mp_set_binary_from_int(mp_number *A, int B);
74 static void mp_set_binary_from_boolean(mp_number *A, int B);
75 static void mp_set_binary_from_scaled(mp_number *A, int B);
76 static void mp_set_binary_from_addition(mp_number *A, mp_number B, mp_number C);
77 static void mp_set_binary_from_substraction (mp_number *A, mp_number B, mp_number C);
78 static void mp_set_binary_from_div(mp_number *A, mp_number B, mp_number C);
79 static void mp_set_binary_from_mul(mp_number *A, mp_number B, mp_number C);
80 static void mp_set_binary_from_int_div(mp_number *A, mp_number B, int C);
81 static void mp_set_binary_from_int_mul(mp_number *A, mp_number B, int C);
82 static void mp_set_binary_from_of_the_way(MP mp, mp_number *A, mp_number t, mp_number B, mp_number C);
83 static void mp_number_negate(mp_number *A);
84 static void mp_number_add(mp_number *A, mp_number B);
85 static void mp_number_substract(mp_number *A, mp_number B);
86 static void mp_number_half(mp_number *A);
87 static void mp_number_halfp(mp_number *A);
88 static void mp_number_double(mp_number *A);
89 static void mp_number_add_scaled(mp_number *A, int B); /* also for negative B */
90 static void mp_number_multiply_int(mp_number *A, int B);
91 static void mp_number_divide_int(mp_number *A, int B);
92 static void mp_binary_abs(mp_number *A);
93 static void mp_number_clone(mp_number *A, mp_number B);
94 static void mp_number_swap(mp_number *A, mp_number *B);
95 static int mp_round_unscaled(mp_number x_orig);
96 static int mp_number_to_int(mp_number A);
97 static int mp_number_to_scaled(mp_number A);
98 static int mp_number_to_boolean(mp_number A);
99 static double mp_number_to_double(mp_number A);
100 static int mp_number_odd(mp_number A);
101 static int mp_number_equal(mp_number A, mp_number B);
102 static int mp_number_greater(mp_number A, mp_number B);
103 static int mp_number_less(mp_number A, mp_number B);
104 static int mp_number_nonequalabs(mp_number A, mp_number B);
105 static void mp_number_floor (mp_number *i);
106 static void mp_binary_fraction_to_round_scaled (mp_number *x);
107 static void mp_binary_number_make_scaled (MP mp, mp_number *r, mp_number p, mp_number q);
108 static void mp_binary_number_make_fraction (MP mp, mp_number *r, mp_number p, mp_number q);
109 static void mp_binary_number_take_fraction (MP mp, mp_number *r, mp_number p, mp_number q);
110 static void mp_binary_number_take_scaled (MP mp, mp_number *r, mp_number p, mp_number q);
111 static void mp_new_number (MP mp, mp_number *n, mp_number_type t) ;
112 static void mp_free_number (MP mp, mp_number *n) ;
113 static void mp_set_binary_from_double(mp_number *A, double B);
114 static void mp_free_binary_math (MP mp);
115 static void mp_binary_set_precision (MP mp);
116 static void mp_check_mpfr_t (MP mp, mpfr_t dec);
117 static int binary_number_check (mpfr_t dec);
118 static char * mp_binnumber_tostring (mpfr_t n);
119 static void init_binary_constants (void);
120 static void free_binary_constants (void);
121 static mpfr_prec_t precision_digits_to_bits(double i);
122 static double precision_bits_to_digits (mpfr_prec_t i);
123
124 @ We do not want special numbers as return values for functions, so:
125
126 @d mpfr_negative_p(a) (mpfr_sgn((a))<0)
127 @d mpfr_positive_p(a) (mpfr_sgn((a))>0)
128 @d checkZero(dec) if (mpfr_zero_p(dec) && mpfr_negative_p(dec)) {
129 mpfr_set_zero(dec,1);
130 }
131
132 @c
133 int binary_number_check (mpfr_t dec)
134 {
135 int test = false;
136 if (!mpfr_number_p(dec)) {
137 test = true;
138 if (mpfr_inf_p(dec)) {
139 mpfr_set(dec, EL_GORDO_mpfr_t, ROUNDING);
140 if (mpfr_negative_p(dec)) {
141 mpfr_neg(dec, dec, ROUNDING);
142 }
143 } else { // Nan
144 mpfr_set_zero(dec,1); /* 1 == positive */
145 }
146 }
147 checkZero(dec);
148 return test;
149 }
150 void mp_check_mpfr_t (MP mp, mpfr_t dec)
151 {
152 mp->arith_error = binary_number_check (dec);
153 }
154
155
156
157
158 @ Precision IO uses |double| because |MPFR_PREC_MAX| overflows int.
159
160 @c
161 static double precision_bits;
162 mpfr_prec_t precision_digits_to_bits (double i)
163 {
164 return i/log10(2);
165 }
166 double precision_bits_to_digits (mpfr_prec_t d)
167 {
168 return d*log10(2);
169 }
170
171
172 @ And these are the ones that {\it are} used elsewhere
173
174 @<Internal library declarations@>=
175 void * mp_initialize_binary_math (MP mp);
176
177 @
178
179 @d unity 1
180 @d two 2
181 @d three 3
182 @d four 4
183 @d half_unit 0.5
184 @d three_quarter_unit 0.75
185 @d coef_bound ((7.0/3.0)*fraction_multiplier) /* |fraction| approximation to 7/3 */
186 @d fraction_threshold 0.04096 /* a |fraction| coefficient less than this is zeroed */
187 @d half_fraction_threshold (fraction_threshold/2) /* half of |fraction_threshold| */
188 @d scaled_threshold 0.000122 /* a |scaled| coefficient less than this is zeroed */
189 @d half_scaled_threshold (scaled_threshold/2) /* half of |scaled_threshold| */
190 @d near_zero_angle (0.0256*angle_multiplier) /* an angle of about 0.0256 */
191 @d p_over_v_threshold 0x80000 /* TODO */
192 @d equation_threshold 0.001
193 @d tfm_warn_threshold 0.0625
194 @d warning_limit pow(2.0,52.0) /* this is a large value that can just be expressed without loss of precision */
195 @d epsilon "1E-52"
196 @d epsilonf pow(2.0,-52.0)
197 @d EL_GORDO "1E1000000" /* the largest value that \MP\ likes. */
198 @d one_third_EL_GORDO (EL_GORDO/3.0)
199
200 @<Declarations@>=
201 static mpfr_t zero;
202 static mpfr_t one;
203 static mpfr_t minusone;
204 static mpfr_t two_mpfr_t;
205 static mpfr_t three_mpfr_t;
206 static mpfr_t four_mpfr_t;
207 static mpfr_t fraction_multiplier_mpfr_t;
208 static mpfr_t angle_multiplier_mpfr_t;
209 static mpfr_t fraction_one_mpfr_t;
210 static mpfr_t fraction_one_plus_mpfr_t;
211 static mpfr_t PI_mpfr_t;
212 static mpfr_t epsilon_mpfr_t;
213 static mpfr_t EL_GORDO_mpfr_t;
214
215 @ @c
216 void init_binary_constants (void) {
217 mpfr_inits2 (precision_bits, one, minusone, zero, two_mpfr_t, three_mpfr_t, four_mpfr_t, fraction_multiplier_mpfr_t,
218 fraction_one_mpfr_t, fraction_one_plus_mpfr_t, angle_multiplier_mpfr_t, PI_mpfr_t,
219 epsilon_mpfr_t, EL_GORDO_mpfr_t, (mpfr_ptr) 0);
220 mpfr_set_si (one, 1, ROUNDING);
221 mpfr_set_si (minusone, -1, ROUNDING);
222 mpfr_set_si (zero, 0, ROUNDING);
223 mpfr_set_si (two_mpfr_t, two, ROUNDING);
224 mpfr_set_si (three_mpfr_t, three, ROUNDING);
225 mpfr_set_si (four_mpfr_t, four, ROUNDING);
226 mpfr_set_si (fraction_multiplier_mpfr_t, fraction_multiplier, ROUNDING);
227 mpfr_set_si (fraction_one_mpfr_t, fraction_one, ROUNDING);
228 mpfr_set_si (fraction_one_plus_mpfr_t, (fraction_one+1), ROUNDING);
229 mpfr_set_si (angle_multiplier_mpfr_t, angle_multiplier, ROUNDING);
230 mpfr_set_str (PI_mpfr_t, PI_STRING, 10, ROUNDING);
231 mpfr_set_str (epsilon_mpfr_t, epsilon, 10, ROUNDING);
232 mpfr_set_str (EL_GORDO_mpfr_t, EL_GORDO, 10, ROUNDING);
233 }
234 void free_binary_constants (void) {
235 mpfr_clears (one, minusone, zero, two_mpfr_t, three_mpfr_t, four_mpfr_t, fraction_multiplier_mpfr_t,
236 fraction_one_mpfr_t, fraction_one_plus_mpfr_t, angle_multiplier_mpfr_t, PI_mpfr_t,
237 epsilon_mpfr_t, EL_GORDO_mpfr_t, (mpfr_ptr) 0);
238 mpfr_free_cache ();
239 }
240
241 @ |precision_max| is limited to 1000, because the precision of already initialized
242 |mpfr_t| numbers cannot be raised, only lowered. The value of 1000.0 is a tradeoff
243 between precision and allocation size / processing speed.
244
245 @d MAX_PRECISION 1000.0
246 @d DEF_PRECISION 34.0
247
248 @c
249 void * mp_initialize_binary_math (MP mp) {
250 math_data *math = (math_data *)mp_xmalloc(mp,1,sizeof(math_data));
251 precision_bits = precision_digits_to_bits(MAX_PRECISION);
252 init_binary_constants();
253 /* alloc */
254 math->allocate = mp_new_number;
255 math->free = mp_free_number;
256 mp_new_number (mp, &math->precision_default, mp_scaled_type);
257 mpfr_set_d(math->precision_default.data.num, DEF_PRECISION, ROUNDING);
258 mp_new_number (mp, &math->precision_max, mp_scaled_type);
259 mpfr_set_d(math->precision_max.data.num, MAX_PRECISION, ROUNDING);
260 mp_new_number (mp, &math->precision_min, mp_scaled_type);
261 /* really should be |precision_bits_to_digits(MPFR_PREC_MIN)| but that produces a horrible number */
262 mpfr_set_d(math->precision_min.data.num, 1.0 , ROUNDING);
263 /* here are the constants for |scaled| objects */
264 mp_new_number (mp, &math->epsilon_t, mp_scaled_type);
265 mpfr_set (math->epsilon_t.data.num, epsilon_mpfr_t, ROUNDING);
266 mp_new_number (mp, &math->inf_t, mp_scaled_type);
267 mpfr_set (math->inf_t.data.num, EL_GORDO_mpfr_t, ROUNDING);
268 mp_new_number (mp, &math->warning_limit_t, mp_scaled_type);
269 mpfr_set_d (math->warning_limit_t.data.num, warning_limit, ROUNDING);
270 mp_new_number (mp, &math->one_third_inf_t, mp_scaled_type);
271 mpfr_div (math->one_third_inf_t.data.num, math->inf_t.data.num, three_mpfr_t, ROUNDING);
272 mp_new_number (mp, &math->unity_t, mp_scaled_type);
273 mpfr_set (math->unity_t.data.num, one, ROUNDING);
274 mp_new_number (mp, &math->two_t, mp_scaled_type);
275 mpfr_set_si(math->two_t.data.num, two, ROUNDING);
276 mp_new_number (mp, &math->three_t, mp_scaled_type);
277 mpfr_set_si(math->three_t.data.num, three, ROUNDING);
278 mp_new_number (mp, &math->half_unit_t, mp_scaled_type);
279 mpfr_set_d(math->half_unit_t.data.num, half_unit, ROUNDING);
280 mp_new_number (mp, &math->three_quarter_unit_t, mp_scaled_type);
281 mpfr_set_d (math->three_quarter_unit_t.data.num, three_quarter_unit, ROUNDING);
282 mp_new_number (mp, &math->zero_t, mp_scaled_type);
283 mpfr_set_zero (math->zero_t.data.num, 1);
284 /* |fractions| */
285 mp_new_number (mp, &math->arc_tol_k, mp_fraction_type);
286 {
287 mpfr_div_si (math->arc_tol_k.data.num, one, 4096, ROUNDING);
288 /* quit when change in arc length estimate reaches this */
289 }
290 mp_new_number (mp, &math->fraction_one_t, mp_fraction_type);
291 mpfr_set_si(math->fraction_one_t.data.num, fraction_one, ROUNDING);
292 mp_new_number (mp, &math->fraction_half_t, mp_fraction_type);
293 mpfr_set_si(math->fraction_half_t.data.num, fraction_half, ROUNDING);
294 mp_new_number (mp, &math->fraction_three_t, mp_fraction_type);
295 mpfr_set_si(math->fraction_three_t.data.num, fraction_three, ROUNDING);
296 mp_new_number (mp, &math->fraction_four_t, mp_fraction_type);
297 mpfr_set_si(math->fraction_four_t.data.num, fraction_four, ROUNDING);
298 /* |angles| */
299 mp_new_number (mp, &math->three_sixty_deg_t, mp_angle_type);
300 mpfr_set_si(math->three_sixty_deg_t.data.num, 360 * angle_multiplier, ROUNDING);
301 mp_new_number (mp, &math->one_eighty_deg_t, mp_angle_type);
302 mpfr_set_si(math->one_eighty_deg_t.data.num, 180 * angle_multiplier, ROUNDING);
303 /* various approximations */
304 mp_new_number (mp, &math->one_k, mp_scaled_type);
305 mpfr_set_d(math->one_k.data.num, 1.0/64, ROUNDING);
306 mp_new_number (mp, &math->sqrt_8_e_k, mp_scaled_type);
307 {
308 mpfr_set_d(math->sqrt_8_e_k.data.num, 112428.82793 / 65536.0, ROUNDING);
309 /* $2^{16}\sqrt{8/e}\approx 112428.82793$ */
310 }
311 mp_new_number (mp, &math->twelve_ln_2_k, mp_fraction_type);
312 {
313 mpfr_set_d(math->twelve_ln_2_k.data.num, 139548959.6165 / 65536.0, ROUNDING);
314 /* $2^{24}\cdot12\ln2\approx139548959.6165$ */
315 }
316 mp_new_number (mp, &math->coef_bound_k, mp_fraction_type);
317 mpfr_set_d(math->coef_bound_k.data.num,coef_bound, ROUNDING);
318 mp_new_number (mp, &math->coef_bound_minus_1, mp_fraction_type);
319 mpfr_set_d(math->coef_bound_minus_1.data.num,coef_bound - 1 / 65536.0, ROUNDING);
320 mp_new_number (mp, &math->twelvebits_3, mp_scaled_type);
321 {
322 mpfr_set_d(math->twelvebits_3.data.num, 1365 / 65536.0, ROUNDING);
323 /* $1365\approx 2^{12}/3$ */
324 }
325 mp_new_number (mp, &math->twentysixbits_sqrt2_t, mp_fraction_type);
326 {
327 mpfr_set_d(math->twentysixbits_sqrt2_t.data.num, 94906265.62 / 65536.0, ROUNDING);
328 /* $2^{26}\sqrt2\approx94906265.62$ */
329 }
330 mp_new_number (mp, &math->twentyeightbits_d_t, mp_fraction_type);
331 {
332 mpfr_set_d(math->twentyeightbits_d_t.data.num, 35596754.69 / 65536.0, ROUNDING);
333 /* $2^{28}d\approx35596754.69$ */
334 }
335 mp_new_number (mp, &math->twentysevenbits_sqrt2_d_t, mp_fraction_type);
336 {
337 mpfr_set_d(math->twentysevenbits_sqrt2_d_t.data.num, 25170706.63 / 65536.0, ROUNDING);
338 /* $2^{27}\sqrt2\,d\approx25170706.63$ */
339 }
340 /* thresholds */
341 mp_new_number (mp, &math->fraction_threshold_t, mp_fraction_type);
342 mpfr_set_d(math->fraction_threshold_t.data.num, fraction_threshold, ROUNDING);
343 mp_new_number (mp, &math->half_fraction_threshold_t, mp_fraction_type);
344 mpfr_set_d(math->half_fraction_threshold_t.data.num, half_fraction_threshold, ROUNDING);
345 mp_new_number (mp, &math->scaled_threshold_t, mp_scaled_type);
346 mpfr_set_d(math->scaled_threshold_t.data.num, scaled_threshold, ROUNDING);
347 mp_new_number (mp, &math->half_scaled_threshold_t, mp_scaled_type);
348 mpfr_set_d(math->half_scaled_threshold_t.data.num, half_scaled_threshold, ROUNDING);
349 mp_new_number (mp, &math->near_zero_angle_t, mp_angle_type);
350 mpfr_set_d(math->near_zero_angle_t.data.num, near_zero_angle, ROUNDING);
351 mp_new_number (mp, &math->p_over_v_threshold_t, mp_fraction_type);
352 mpfr_set_d(math->p_over_v_threshold_t.data.num, p_over_v_threshold, ROUNDING);
353 mp_new_number (mp, &math->equation_threshold_t, mp_scaled_type);
354 mpfr_set_d(math->equation_threshold_t.data.num, equation_threshold, ROUNDING);
355 mp_new_number (mp, &math->tfm_warn_threshold_t, mp_scaled_type);
356 mpfr_set_d(math->tfm_warn_threshold_t.data.num, tfm_warn_threshold, ROUNDING);
357 /* functions */
358 math->from_int = mp_set_binary_from_int;
359 math->from_boolean = mp_set_binary_from_boolean;
360 math->from_scaled = mp_set_binary_from_scaled;
361 math->from_double = mp_set_binary_from_double;
362 math->from_addition = mp_set_binary_from_addition;
363 math->from_substraction = mp_set_binary_from_substraction;
364 math->from_oftheway = mp_set_binary_from_of_the_way;
365 math->from_div = mp_set_binary_from_div;
366 math->from_mul = mp_set_binary_from_mul;
367 math->from_int_div = mp_set_binary_from_int_div;
368 math->from_int_mul = mp_set_binary_from_int_mul;
369 math->negate = mp_number_negate;
370 math->add = mp_number_add;
371 math->substract = mp_number_substract;
372 math->half = mp_number_half;
373 math->halfp = mp_number_halfp;
374 math->do_double = mp_number_double;
375 math->abs = mp_binary_abs;
376 math->clone = mp_number_clone;
377 math->swap = mp_number_swap;
378 math->add_scaled = mp_number_add_scaled;
379 math->multiply_int = mp_number_multiply_int;
380 math->divide_int = mp_number_divide_int;
381 math->to_boolean = mp_number_to_boolean;
382 math->to_scaled = mp_number_to_scaled;
383 math->to_double = mp_number_to_double;
384 math->to_int = mp_number_to_int;
385 math->odd = mp_number_odd;
386 math->equal = mp_number_equal;
387 math->less = mp_number_less;
388 math->greater = mp_number_greater;
389 math->nonequalabs = mp_number_nonequalabs;
390 math->round_unscaled = mp_round_unscaled;
391 math->floor_scaled = mp_number_floor;
392 math->fraction_to_round_scaled = mp_binary_fraction_to_round_scaled;
393 math->make_scaled = mp_binary_number_make_scaled;
394 math->make_fraction = mp_binary_number_make_fraction;
395 math->take_fraction = mp_binary_number_take_fraction;
396 math->take_scaled = mp_binary_number_take_scaled;
397 math->velocity = mp_binary_velocity;
398 math->n_arg = mp_binary_n_arg;
399 math->m_log = mp_binary_m_log;
400 math->m_exp = mp_binary_m_exp;
401 math->m_unif_rand = mp_binary_m_unif_rand;
402 math->m_norm_rand = mp_binary_m_norm_rand;
403 math->pyth_add = mp_binary_pyth_add;
404 math->pyth_sub = mp_binary_pyth_sub;
405 math->fraction_to_scaled = mp_number_fraction_to_scaled;
406 math->scaled_to_fraction = mp_number_scaled_to_fraction;
407 math->scaled_to_angle = mp_number_scaled_to_angle;
408 math->angle_to_scaled = mp_number_angle_to_scaled;
409 math->init_randoms = mp_init_randoms;
410 math->sin_cos = mp_binary_sin_cos;
411 math->slow_add = mp_binary_slow_add;
412 math->sqrt = mp_binary_square_rt;
413 math->print = mp_binary_print_number;
414 math->tostring = mp_binary_number_tostring;
415 math->modulo = mp_binary_number_modulo;
416 math->ab_vs_cd = mp_ab_vs_cd;
417 math->crossing_point = mp_binary_crossing_point;
418 math->scan_numeric = mp_binary_scan_numeric_token;
419 math->scan_fractional = mp_binary_scan_fractional_token;
420 math->free_math = mp_free_binary_math;
421 math->set_precision = mp_binary_set_precision;
422 return (void *)math;
423 }
424
425 void mp_binary_set_precision (MP mp) {
426 double d = mpfr_get_d(internal_value (mp_number_precision).data.num, ROUNDING);
427 precision_bits = precision_digits_to_bits(d);
428 }
429
430 void mp_free_binary_math (MP mp) {
431 free_number (((math_data *)mp->math)->three_sixty_deg_t);
432 free_number (((math_data *)mp->math)->one_eighty_deg_t);
433 free_number (((math_data *)mp->math)->fraction_one_t);
434 free_number (((math_data *)mp->math)->zero_t);
435 free_number (((math_data *)mp->math)->half_unit_t);
436 free_number (((math_data *)mp->math)->three_quarter_unit_t);
437 free_number (((math_data *)mp->math)->unity_t);
438 free_number (((math_data *)mp->math)->two_t);
439 free_number (((math_data *)mp->math)->three_t);
440 free_number (((math_data *)mp->math)->one_third_inf_t);
441 free_number (((math_data *)mp->math)->inf_t);
442 free_number (((math_data *)mp->math)->warning_limit_t);
443 free_number (((math_data *)mp->math)->one_k);
444 free_number (((math_data *)mp->math)->sqrt_8_e_k);
445 free_number (((math_data *)mp->math)->twelve_ln_2_k);
446 free_number (((math_data *)mp->math)->coef_bound_k);
447 free_number (((math_data *)mp->math)->coef_bound_minus_1);
448 free_number (((math_data *)mp->math)->fraction_threshold_t);
449 free_number (((math_data *)mp->math)->half_fraction_threshold_t);
450 free_number (((math_data *)mp->math)->scaled_threshold_t);
451 free_number (((math_data *)mp->math)->half_scaled_threshold_t);
452 free_number (((math_data *)mp->math)->near_zero_angle_t);
453 free_number (((math_data *)mp->math)->p_over_v_threshold_t);
454 free_number (((math_data *)mp->math)->equation_threshold_t);
455 free_number (((math_data *)mp->math)->tfm_warn_threshold_t);
456 free_binary_constants();
457 free(mp->math);
458 }
459
460 @ Creating an destroying |mp_number| objects
461
462 @ @c
463 void mp_new_number (MP mp, mp_number *n, mp_number_type t) {
464 (void)mp;
465 n->data.num = mp_xmalloc(mp,1,sizeof(mpfr_t));
466 mpfr_init2 ((mpfr_ptr)(n->data.num), precision_bits);
467 mpfr_set_zero((mpfr_ptr)(n->data.num),1); /* 1 == positive */
468 n->type = t;
469 }
470
471 @
472
473 @c
474 void mp_free_number (MP mp, mp_number *n) {
475 (void)mp;
476 if (n->data.num) {
477 mpfr_clear (n->data.num);
478 n->data.num = NULL;
479 }
480 n->type = mp_nan_type;
481 }
482
483 @ Here are the low-level functions on |mp_number| items, setters first.
484
485 @c
486 void mp_set_binary_from_int(mp_number *A, int B) {
487 mpfr_set_si(A->data.num,B, ROUNDING);
488 }
489 void mp_set_binary_from_boolean(mp_number *A, int B) {
490 mpfr_set_si(A->data.num,B, ROUNDING);
491 }
492 void mp_set_binary_from_scaled(mp_number *A, int B) {
493 mpfr_set_si(A->data.num, B, ROUNDING);
494 mpfr_div_si(A->data.num, A->data.num, 65536, ROUNDING);
495 }
496 void mp_set_binary_from_double(mp_number *A, double B) {
497 mpfr_set_d(A->data.num, B, ROUNDING);
498 }
499 void mp_set_binary_from_addition(mp_number *A, mp_number B, mp_number C) {
500 mpfr_add(A->data.num,B.data.num,C.data.num, ROUNDING);
501 }
502 void mp_set_binary_from_substraction (mp_number *A, mp_number B, mp_number C) {
503 mpfr_sub(A->data.num,B.data.num,C.data.num, ROUNDING);
504 }
505 void mp_set_binary_from_div(mp_number *A, mp_number B, mp_number C) {
506 mpfr_div(A->data.num,B.data.num,C.data.num, ROUNDING);
507 }
508 void mp_set_binary_from_mul(mp_number *A, mp_number B, mp_number C) {
509 mpfr_mul(A->data.num,B.data.num,C.data.num, ROUNDING);
510 }
511 void mp_set_binary_from_int_div(mp_number *A, mp_number B, int C) {
512 mpfr_div_si(A->data.num,B.data.num,C, ROUNDING);
513 }
514 void mp_set_binary_from_int_mul(mp_number *A, mp_number B, int C) {
515 mpfr_mul_si(A->data.num,B.data.num, C, ROUNDING);
516 }
517 void mp_set_binary_from_of_the_way(MP mp, mp_number *A, mp_number t, mp_number B, mp_number C) {
518 mpfr_t c, r1;
519 mpfr_init2(c, precision_bits);
520 mpfr_init2(r1, precision_bits);
521 mpfr_sub (c,B.data.num, C.data.num, ROUNDING);
522 mp_binary_take_fraction(mp, r1, c, t.data.num);
523 mpfr_sub (A->data.num, B.data.num, r1, ROUNDING);
524 mpfr_clear(c);
525 mpfr_clear(r1);
526 mp_check_mpfr_t(mp, A->data.num);
527 }
528 void mp_number_negate(mp_number *A) {
529 mpfr_neg (A->data.num, A->data.num, ROUNDING);
530 checkZero((mpfr_ptr)A->data.num);
531 }
532 void mp_number_add(mp_number *A, mp_number B) {
533 mpfr_add (A->data.num,A->data.num,B.data.num, ROUNDING);
534 }
535 void mp_number_substract(mp_number *A, mp_number B) {
536 mpfr_sub (A->data.num,A->data.num,B.data.num, ROUNDING);
537 }
538 void mp_number_half(mp_number *A) {
539 mpfr_div_si(A->data.num, A->data.num, 2, ROUNDING);
540 }
541 void mp_number_halfp(mp_number *A) {
542 mpfr_div_si(A->data.num,A->data.num, 2, ROUNDING);
543 }
544 void mp_number_double(mp_number *A) {
545 mpfr_mul_si(A->data.num,A->data.num, 2, ROUNDING);
546 }
547 void mp_number_add_scaled(mp_number *A, int B) { /* also for negative B */
548 mpfr_add_d (A->data.num,A->data.num, B/65536.0, ROUNDING);
549 }
550 void mp_number_multiply_int(mp_number *A, int B) {
551 mpfr_mul_si(A->data.num,A->data.num, B, ROUNDING);
552 }
553 void mp_number_divide_int(mp_number *A, int B) {
554 mpfr_div_si(A->data.num,A->data.num, B, ROUNDING);
555 }
556 void mp_binary_abs(mp_number *A) {
557 mpfr_abs(A->data.num, A->data.num, ROUNDING);
558 }
559 void mp_number_clone(mp_number *A, mp_number B) {
560 mpfr_prec_round (A->data.num, precision_bits, ROUNDING);
561 mpfr_set(A->data.num, (mpfr_ptr)B.data.num, ROUNDING);
562 }
563 void mp_number_swap(mp_number *A, mp_number *B) {
564 mpfr_swap(A->data.num, B->data.num);
565 }
566 void mp_number_fraction_to_scaled (mp_number *A) {
567 A->type = mp_scaled_type;
568 mpfr_div (A->data.num, A->data.num, fraction_multiplier_mpfr_t, ROUNDING);
569 }
570 void mp_number_angle_to_scaled (mp_number *A) {
571 A->type = mp_scaled_type;
572 mpfr_div (A->data.num, A->data.num, angle_multiplier_mpfr_t, ROUNDING);
573 }
574 void mp_number_scaled_to_fraction (mp_number *A) {
575 A->type = mp_fraction_type;
576 mpfr_mul (A->data.num, A->data.num, fraction_multiplier_mpfr_t, ROUNDING);
577 }
578 void mp_number_scaled_to_angle (mp_number *A) {
579 A->type = mp_angle_type;
580 mpfr_mul(A->data.num, A->data.num, angle_multiplier_mpfr_t, ROUNDING);
581 }
582
583
584 @* Query functions
585
586 @ Convert a number to a scaled value. |decNumberToInt32| is not
587 able to make this conversion properly, so instead we are using
588 |decNumberToDouble| and a typecast. Bad!
589
590 @c
591 int mp_number_to_scaled(mp_number A) {
592 double v = mpfr_get_d (A.data.num, ROUNDING);
593 return (int)(v * 65536.0);
594 }
595
596 @
597
598 @d odd(A) (abs(A)%2==1)
599
600 @c
601 int mp_number_to_int(mp_number A) {
602 int32_t result = 0;
603 if (mpfr_fits_sint_p(A.data.num, ROUNDING)) {
604 result = mpfr_get_si(A.data.num, ROUNDING);
605 }
606 return result;
607 }
608 int mp_number_to_boolean(mp_number A) {
609 int32_t result = 0;
610 if (mpfr_fits_sint_p(A.data.num, ROUNDING)) {
611 result = mpfr_get_si(A.data.num, ROUNDING);
612 }
613 return result;
614 }
615 double mp_number_to_double(mp_number A) {
616 double res = 0.0;
617 if (mpfr_number_p (A.data.num)) {
618 res = mpfr_get_d(A.data.num, ROUNDING);
619 }
620 return res;
621 }
622 int mp_number_odd(mp_number A) {
623 return odd(mp_number_to_int(A));
624 }
625 int mp_number_equal(mp_number A, mp_number B) {
626 return mpfr_equal_p(A.data.num,B.data.num);
627 }
628 int mp_number_greater(mp_number A, mp_number B) {
629 return mpfr_greater_p(A.data.num,B.data.num);
630 }
631 int mp_number_less(mp_number A, mp_number B) {
632 return mpfr_less_p(A.data.num,B.data.num);
633 }
634 int mp_number_nonequalabs(mp_number A, mp_number B) {
635 return !(mpfr_cmpabs(A.data.num, B.data.num)==0);
636 }
637
638 @ Fixed-point arithmetic is done on {\sl scaled integers\/} that are multiples
639 of $2^{-16}$. In other words, a binary point is assumed to be sixteen bit
640 positions from the right end of a binary computer word.
641
642 @ One of \MP's most common operations is the calculation of
643 $\lfloor{a+b\over2}\rfloor$,
644 the midpoint of two given integers |a| and~|b|. The most decent way to do
645 this is to write `|(a+b)/2|'; but on many machines it is more efficient
646 to calculate `|(a+b)>>1|'.
647
648 Therefore the midpoint operation will always be denoted by `|half(a+b)|'
649 in this program. If \MP\ is being implemented with languages that permit
650 binary shifting, the |half| macro should be changed to make this operation
651 as efficient as possible. Since some systems have shift operators that can
652 only be trusted to work on positive numbers, there is also a macro |halfp|
653 that is used only when the quantity being halved is known to be positive
654 or zero.
655
656 @ Here is a procedure analogous to |print_int|. The current version
657 is fairly stupid, and it is not round-trip safe, but this is good
658 enough for a beta test.
659
660 @c
661 char * mp_binnumber_tostring (mpfr_t n) {
662 char *str = NULL, *buffer = NULL;
663 mpfr_exp_t exp = 0;
664 int neg = 0;
665 if ((str = mpfr_get_str (NULL, &exp, 10, 0, n, ROUNDING))>0) {
666 int numprecdigits = precision_bits_to_digits(precision_bits);
667 if (*str == '-') {
668 neg = 1;
669 }
670 while (strlen(str)>0 && *(str+strlen(str)-1) == '0' ) {
671 *(str+strlen(str)-1) = '\0'; /* get rid of trailing zeroes */
672 }
673 buffer = malloc(strlen(str)+13+numprecdigits+1);
674 /* the buffer should also fit at least strlen("E+%d", exp) or (numprecdigits-2) worth of zeroes,
675 * because with numprecdigits == 33, the str for "1E32" will be "1", and needing 32 extra zeroes,
676 * and the decimal dot. To avoid miscalculations by myself, it is safer to add these
677 * three together.
678 */
679 if (buffer) {
680 int i = 0, j = 0;
681 if (neg) {
682 buffer[i++] = '-';
683 j = 1;
684 }
685 if (strlen(str+j) == 0) {
686 buffer[i++] = '0';
687 } else {
688 /* non-zero */
689 if (exp<=numprecdigits && exp > -6) {
690 if (exp>0) {
691 buffer[i++] = str[j++];
692 while (--exp>0) {
693 buffer[i++] = (str[j] ? str[j++] : '0');
694 }
695 if (str[j]) {
696 buffer[i++] = '.';
697 while (str[j]) {
698 buffer[i++] = str[j++];
699 }
700 }
701 } else {
702 int absexp;
703 buffer[i++] = '0';
704 buffer[i++] = '.';
705 absexp = -exp;
706 while (absexp-- > 0) {
707 buffer[i++] = '0';
708 }
709 while (str[j]) {
710 buffer[i++] = str[j++];
711 }
712 }
713 } else {
714 buffer[i++] = str[j++];
715 if (str[j]) {
716 buffer[i++] = '.';
717 while (str[j]) {
718 buffer[i++] = str[j++];
719 }
720 }
721 {
722 char msg[256];
723 int k = 0;
724 mp_snprintf (msg, 256, "%s%d", (exp>0?"+":""), (int)(exp>0 ? (exp-1) : (exp-1)));
725 buffer[i++] = 'E';
726 while (msg[k]) {
727 buffer[i++] = msg[k++];
728 }
729 }
730 }
731 }
732 buffer[i++] = '\0';
733 }
734 mpfr_free_str(str);
735 }
736 return buffer;
737 }
738 char * mp_binary_number_tostring (MP mp, mp_number n) {
739 return mp_binnumber_tostring(n.data.num);
740 }
741
742
743 @ @c
744 void mp_binary_print_number (MP mp, mp_number n) {
745 char *str = mp_binary_number_tostring(mp, n);
746 mp_print (mp, str);
747 free (str);
748 }
749
750
751
752
753 @ Addition is not always checked to make sure that it doesn't overflow,
754 but in places where overflow isn't too unlikely the |slow_add| routine
755 is used.
756
757 @c
758 void mp_binary_slow_add (MP mp, mp_number *ret, mp_number A, mp_number B) {
759 mpfr_add(ret->data.num,A.data.num,B.data.num, ROUNDING);
760 }
761
762 @ The |make_fraction| routine produces the |fraction| equivalent of
763 |p/q|, given integers |p| and~|q|; it computes the integer
764 $f=\lfloor2^{28}p/q+{1\over2}\rfloor$, when $p$ and $q$ are
765 positive. If |p| and |q| are both of the same scaled type |t|,
766 the ``type relation'' |make_fraction(t,t)=fraction| is valid;
767 and it's also possible to use the subroutine ``backwards,'' using
768 the relation |make_fraction(t,fraction)=t| between scaled types.
769
770 If the result would have magnitude $2^{31}$ or more, |make_fraction|
771 sets |arith_error:=true|. Most of \MP's internal computations have
772 been designed to avoid this sort of error.
773
774 If this subroutine were programmed in assembly language on a typical
775 machine, we could simply compute |(@t$2^{28}$@>*p)div q|, since a
776 double-precision product can often be input to a fixed-point division
777 instruction. But when we are restricted to int-eger arithmetic it
778 is necessary either to resort to multiple-precision maneuvering
779 or to use a simple but slow iteration. The multiple-precision technique
780 would be about three times faster than the code adopted here, but it
781 would be comparatively long and tricky, involving about sixteen
782 additional multiplications and divisions.
783
784 This operation is part of \MP's ``inner loop''; indeed, it will
785 consume nearly 10\pct! of the running time (exclusive of input and output)
786 if the code below is left unchanged. A machine-dependent recoding
787 will therefore make \MP\ run faster. The present implementation
788 is highly portable, but slow; it avoids multiplication and division
789 except in the initial stage. System wizards should be careful to
790 replace it with a routine that is guaranteed to produce identical
791 results in all cases.
792 @^system dependencies@>
793
794 As noted below, a few more routines should also be replaced by machine-dependent
795 code, for efficiency. But when a procedure is not part of the ``inner loop,''
796 such changes aren't advisable; simplicity and robustness are
797 preferable to trickery, unless the cost is too high.
798 @^inner loop@>
799
800 @c
801 void mp_binary_make_fraction (MP mp, mpfr_t ret, mpfr_t p, mpfr_t q) {
802 mpfr_div (ret, p, q, ROUNDING);
803 mp_check_mpfr_t(mp, ret);
804 mpfr_mul (ret, ret, fraction_multiplier_mpfr_t, ROUNDING);
805 }
806 void mp_binary_number_make_fraction (MP mp, mp_number *ret, mp_number p, mp_number q) {
807 mp_binary_make_fraction (mp, ret->data.num, p.data.num, q.data.num);
808 }
809
810 @ @<Declarations@>=
811 void mp_binary_make_fraction (MP mp, mpfr_t ret, mpfr_t p, mpfr_t q);
812
813 @ The dual of |make_fraction| is |take_fraction|, which multiplies a
814 given integer~|q| by a fraction~|f|. When the operands are positive, it
815 computes $p=\lfloor qf/2^{28}+{1\over2}\rfloor$, a symmetric function
816 of |q| and~|f|.
817
818 This routine is even more ``inner loopy'' than |make_fraction|;
819 the present implementation consumes almost 20\pct! of \MP's computation
820 time during typical jobs, so a machine-language substitute is advisable.
821 @^inner loop@> @^system dependencies@>
822
823 @c
824 void mp_binary_take_fraction (MP mp, mpfr_t ret, mpfr_t p, mpfr_t q) {
825 mpfr_mul(ret, p, q, ROUNDING);
826 mpfr_div(ret, ret, fraction_multiplier_mpfr_t, ROUNDING);
827 }
828 void mp_binary_number_take_fraction (MP mp, mp_number *ret, mp_number p, mp_number q) {
829 mp_binary_take_fraction (mp, ret->data.num, p.data.num, q.data.num);
830 }
831
832 @ @<Declarations@>=
833 void mp_binary_take_fraction (MP mp, mpfr_t ret, mpfr_t p, mpfr_t q);
834
835 @ When we want to multiply something by a |scaled| quantity, we use a scheme
836 analogous to |take_fraction| but with a different scaling.
837 Given positive operands, |take_scaled|
838 computes the quantity $p=\lfloor qf/2^{16}+{1\over2}\rfloor$.
839
840 Once again it is a good idea to use a machine-language replacement if
841 possible; otherwise |take_scaled| will use more than 2\pct! of the running time
842 when the Computer Modern fonts are being generated.
843 @^inner loop@>
844
845 @c
846 void mp_binary_number_take_scaled (MP mp, mp_number *ret, mp_number p_orig, mp_number q_orig) {
847 mpfr_mul(ret->data.num, p_orig.data.num, q_orig.data.num, ROUNDING);
848 }
849
850
851 @ For completeness, there's also |make_scaled|, which computes a
852 quotient as a |scaled| number instead of as a |fraction|.
853 In other words, the result is $\lfloor2^{16}p/q+{1\over2}\rfloor$, if the
854 operands are positive. \ (This procedure is not used especially often,
855 so it is not part of \MP's inner loop.)
856
857 @c
858 void mp_binary_number_make_scaled (MP mp, mp_number *ret, mp_number p_orig, mp_number q_orig) {
859 mpfr_div(ret->data.num, p_orig.data.num, q_orig.data.num, ROUNDING);
860 mp_check_mpfr_t(mp, ret->data.num);
861 }
862
863 @
864 @d halfp(A) (integer)((unsigned)(A) >> 1)
865
866 @* Scanning numbers in the input
867
868 The definitions below are temporarily here
869
870 @d set_cur_cmd(A) mp->cur_mod_->type=(A)
871 @d set_cur_mod(A) mpfr_set((mpfr_ptr)(mp->cur_mod_->data.n.data.num),A, ROUNDING)
872
873 @<Declarations...@>=
874 static void mp_wrapup_numeric_token(MP mp, unsigned char *start, unsigned char *stop);
875
876 @ The check of the precision is based on the article "27 Bits are not enough for 8-Digit accuracy"
877 @ by Bennet Goldberg which roughly says that
878 @ given $p$ digits in base 10 and $q$ digits in base 2,
879 @ conversion from base 10 round-trip through base 2 if and only if $10^p < $2^{q-1}$.
880 @ In our case $p/\log_{10}2 + 1 < q$, or $q\geq a$
881 @ where $q$ is the current precision in bits and $a=\left\lceil p/\log_{10}2 + 1\right\rceil$.
882 @ Therefore if $a>q$ the number could be too precise and we emit a warning.
883 @d too_precise(a) (a>precision_bits)
884 @c
885 void mp_wrapup_numeric_token(MP mp, unsigned char *start, unsigned char *stop) {
886 int invalid = 0;
887 mpfr_t result;
888 size_t l = stop-start+1;
889 unsigned long lp, lpbit;
890 /*size_t lp = l; */
891 char *buf = mp_xmalloc(mp, l+1, 1);
892 char *bufp = buf;
893 buf[l] = '\0';
894 mpfr_init2(result, precision_bits);
895 (void)strncpy(buf,(const char *)start, l);
896 invalid = mpfr_set_str(result,buf, 10, ROUNDING);
897 //fprintf(stdout,"scan of [%s] produced %s, ", buf, mp_binnumber_tostring(result));
898 lp = (unsigned long) l;
899 /* strip leading - or + or 0 or .*/
900 if ( (*bufp=='-') || (*bufp=='+') || (*bufp=='0') || (*bufp=='.') ) { lp--; bufp++;}
901 /* strip also . */
902 lp = strchr(bufp,'.') ? lp-1: lp;
903 /* strip also trailing 0s */
904 bufp = buf+l-1;
905 while(*bufp == '0') {bufp--; lp=( ((lp==0)||(lp==1))?1:lp-1);}
906 /* at least one digit, even if the number is 0 */
907 lp = lp>0? lp: 1;
908 /* bits needed for buf */
909 lpbit = (unsigned long)ceil(lp/log10(2)+1);
910 free(buf);
911 bufp = NULL;
912 if (invalid == 0) {
913 set_cur_mod(result);
914 // fprintf(stdout,"mod=%s\n", mp_binary_number_tostring(mp,mp->cur_mod_->data.n));
915 if (too_precise(lpbit)) {
916 if (mpfr_positive_p((mpfr_ptr)(internal_value (mp_warning_check).data.num)) &&
917 (mp->scanner_status != tex_flushing)) {
918 char msg[256];
919 const char *hlp[] = {"Continue and I'll try to cope",
920 "with that value; but it might be dangerous.",
921 "(Set warningcheck:=0 to suppress this message.)",
922 NULL };
923 mp_snprintf (msg, 256, "Number is too precise (%d vs. numberprecision = %f)", (unsigned int)lp,mpfr_get_d(internal_value (mp_number_precision).data.num, ROUNDING));
924 @.Number is too large@>;
925 mp_error (mp, msg, hlp, true);
926 }
927 }
928 } else if (mp->scanner_status != tex_flushing) {
929 const char *hlp[] = {"I could not handle this number specification",
930 "probably because it is out of range. Error:",
931 "",
932 NULL };
933 hlp[2] = strerror(errno);
934 mp_error (mp, "Enormous number has been reduced.", hlp, false);
935 @.Enormous number...@>;
936 set_cur_mod((mpfr_ptr)(((math_data *)(mp->math))->inf_t.data.num));
937 }
938 set_cur_cmd((mp_variable_type)mp_numeric_token);
939 mpfr_clear(result);
940 }
941
942 @ @c
943 static void find_exponent (MP mp) {
944 if (mp->buffer[mp->cur_input.loc_field] == 'e' ||
945 mp->buffer[mp->cur_input.loc_field] == 'E') {
946 mp->cur_input.loc_field++;
947 if (!(mp->buffer[mp->cur_input.loc_field] == '+' ||
948 mp->buffer[mp->cur_input.loc_field] == '-' ||
949 mp->char_class[mp->buffer[mp->cur_input.loc_field]] == digit_class)) {
950 mp->cur_input.loc_field--;
951 return;
952 }
953 if (mp->buffer[mp->cur_input.loc_field] == '+' ||
954 mp->buffer[mp->cur_input.loc_field] == '-') {
955 mp->cur_input.loc_field++;
956 }
957 while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == digit_class) {
958 mp->cur_input.loc_field++;
959 }
960 }
961 }
962 void mp_binary_scan_fractional_token (MP mp, int n) { /* n: scaled */
963 unsigned char *start = &mp->buffer[mp->cur_input.loc_field -1];
964 unsigned char *stop;
965 while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == digit_class) {
966 mp->cur_input.loc_field++;
967 }
968 find_exponent(mp);
969 stop = &mp->buffer[mp->cur_input.loc_field-1];
970 mp_wrapup_numeric_token (mp, start, stop);
971 }
972
973
974 @ We just have to collect bytes.
975
976 @c
977 void mp_binary_scan_numeric_token (MP mp, int n) { /* n: scaled */
978 unsigned char *start = &mp->buffer[mp->cur_input.loc_field -1];
979 unsigned char *stop;
980 while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == digit_class) {
981 mp->cur_input.loc_field++;
982 }
983 if (mp->buffer[mp->cur_input.loc_field] == '.' &&
984 mp->buffer[mp->cur_input.loc_field+1] != '.') {
985 mp->cur_input.loc_field++;
986 while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == digit_class) {
987 mp->cur_input.loc_field++;
988 }
989 }
990 find_exponent(mp);
991 stop = &mp->buffer[mp->cur_input.loc_field-1];
992 mp_wrapup_numeric_token (mp, start, stop);
993 }
994
995 @ The |scaled| quantities in \MP\ programs are generally supposed to be
996 less than $2^{12}$ in absolute value, so \MP\ does much of its internal
997 arithmetic with 28~significant bits of precision. A |fraction| denotes
998 a scaled integer whose binary point is assumed to be 28 bit positions
999 from the right.
1000
1001 @d fraction_half (fraction_multiplier/2)
1002 @d fraction_one (1*fraction_multiplier)
1003 @d fraction_two (2*fraction_multiplier)
1004 @d fraction_three (3*fraction_multiplier)
1005 @d fraction_four (4*fraction_multiplier)
1006
1007 @ Here is a typical example of how the routines above can be used.
1008 It computes the function
1009 $${1\over3\tau}f(\theta,\phi)=
1010 {\tau^{-1}\bigl(2+\sqrt2\,(\sin\theta-{1\over16}\sin\phi)
1011 (\sin\phi-{1\over16}\sin\theta)(\cos\theta-\cos\phi)\bigr)\over
1012 3\,\bigl(1+{1\over2}(\sqrt5-1)\cos\theta+{1\over2}(3-\sqrt5\,)\cos\phi\bigr)},$$
1013 where $\tau$ is a |scaled| ``tension'' parameter. This is \MP's magic
1014 fudge factor for placing the first control point of a curve that starts
1015 at an angle $\theta$ and ends at an angle $\phi$ from the straight path.
1016 (Actually, if the stated quantity exceeds 4, \MP\ reduces it to~4.)
1017
1018 The trigonometric quantity to be multiplied by $\sqrt2$ is less than $\sqrt2$.
1019 (It's a sum of eight terms whose absolute values can be bounded using
1020 relations such as $\sin\theta\cos\theta\L{1\over2}$.) Thus the numerator
1021 is positive; and since the tension $\tau$ is constrained to be at least
1022 $3\over4$, the numerator is less than $16\over3$. The denominator is
1023 nonnegative and at most~6.
1024
1025 The angles $\theta$ and $\phi$ are given implicitly in terms of |fraction|
1026 arguments |st|, |ct|, |sf|, and |cf|, representing $\sin\theta$, $\cos\theta$,
1027 $\sin\phi$, and $\cos\phi$, respectively.
1028
1029 @c
1030 void mp_binary_velocity (MP mp, mp_number *ret, mp_number st, mp_number ct, mp_number sf,
1031 mp_number cf, mp_number t) {
1032 mpfr_t acc, num, denom; /* registers for intermediate calculations */
1033 mpfr_t r1, r2;
1034 mpfr_t arg1, arg2;
1035 mpfr_t i16, fone, fhalf, ftwo, sqrtfive;
1036 mpfr_inits2 (precision_bits, acc, num, denom, r1, r2, arg1, arg2, i16, fone, fhalf, ftwo, sqrtfive, (mpfr_ptr)0);
1037 mpfr_set_si(i16, 16, ROUNDING);
1038 mpfr_set_si(fone, fraction_one, ROUNDING);
1039 mpfr_set_si(fhalf, fraction_half, ROUNDING);
1040 mpfr_set_si(ftwo, fraction_two, ROUNDING);
1041 mpfr_set_si(sqrtfive, 5, ROUNDING);
1042 mpfr_sqrt (sqrtfive, sqrtfive, ROUNDING);
1043 mpfr_div (arg1,sf.data.num, i16, ROUNDING); // arg1 = sf / 16
1044 mpfr_sub (arg1,st.data.num, arg1, ROUNDING); // arg1 = st - arg1
1045 mpfr_div (arg2,st.data.num, i16, ROUNDING); // arg2 = st / 16
1046 mpfr_sub (arg2,sf.data.num, arg2, ROUNDING); // arg2 = sf - arg2
1047 mp_binary_take_fraction (mp, acc, arg1, arg2); // acc = (arg1 * arg2) / fmul
1048
1049 mpfr_set (arg1, acc, ROUNDING);
1050 mpfr_sub (arg2, ct.data.num, cf.data.num, ROUNDING); // arg2 = ct - cf
1051 mp_binary_take_fraction (mp, acc, arg1, arg2); // acc = (arg1 * arg2 ) / fmul
1052
1053 mpfr_sqrt(arg1, two_mpfr_t, ROUNDING); // arg1 = sqrt(2)
1054 mpfr_mul(arg1, arg1, fone, ROUNDING); // arg1 = arg1 * fmul
1055 mp_binary_take_fraction (mp, r1, acc, arg1); // r1 = (acc * arg1) / fmul
1056 mpfr_add(num, ftwo, r1, ROUNDING); // num = ftwo + r1
1057
1058 mpfr_sub(arg1,sqrtfive, one, ROUNDING); // arg1 = sqrt(5) - 1
1059 mpfr_mul(arg1,arg1,fhalf, ROUNDING); // arg1 = arg1 * fmul/2
1060 mpfr_mul(arg1,arg1,three_mpfr_t, ROUNDING); // arg1 = arg1 * 3
1061
1062 mpfr_sub(arg2,three_mpfr_t, sqrtfive, ROUNDING); // arg2 = 3 - sqrt(5)
1063 mpfr_mul(arg2,arg2,fhalf, ROUNDING); // arg2 = arg2 * fmul/2
1064 mpfr_mul(arg2,arg2,three_mpfr_t, ROUNDING); // arg2 = arg2 * 3
1065 mp_binary_take_fraction (mp, r1, ct.data.num, arg1) ; // r1 = (ct * arg1) / fmul
1066 mp_binary_take_fraction (mp, r2, cf.data.num, arg2); // r2 = (cf * arg2) / fmul
1067
1068 mpfr_set_si(denom, fraction_three, ROUNDING); // denom = 3fmul
1069 mpfr_add(denom, denom, r1, ROUNDING); // denom = denom + r1
1070 mpfr_add(denom, denom, r2, ROUNDING); // denom = denom + r2
1071
1072 if (!mpfr_equal_p(t.data.num, one)) { // t != 1
1073 mpfr_div(num, num, t.data.num, ROUNDING); // num = num / t
1074 }
1075 mpfr_set(r2, num, ROUNDING); // r2 = num / 4
1076 mpfr_div(r2, r2, four_mpfr_t, ROUNDING);
1077 if (mpfr_less_p(denom,r2)) { // num/4 >= denom => denom < num/4
1078 mpfr_set_si(ret->data.num,fraction_four, ROUNDING);
1079 } else {
1080 mp_binary_make_fraction (mp, ret->data.num, num, denom);
1081 }
1082 mpfr_clears (acc, num, denom, r1, r2, arg1, arg2, i16, fone, fhalf, ftwo, sqrtfive, (mpfr_ptr)0);
1083 mp_check_mpfr_t(mp, ret->data.num);
1084 }
1085
1086
1087 @ The following somewhat different subroutine tests rigorously if $ab$ is
1088 greater than, equal to, or less than~$cd$,
1089 given integers $(a,b,c,d)$. In most cases a quick decision is reached.
1090 The result is $+1$, 0, or~$-1$ in the three respective cases.
1091
1092 @c
1093 void mp_ab_vs_cd (MP mp, mp_number *ret, mp_number a_orig, mp_number b_orig, mp_number c_orig, mp_number d_orig) {
1094 mpfr_t q, r, test; /* temporary registers */
1095 mpfr_t a, b, c, d;
1096 int cmp = 0;
1097 (void)mp;
1098 mpfr_inits2(precision_bits, q,r,test,a,b,c,d,(mpfr_ptr)0);
1099 mpfr_set(a, (mpfr_ptr)a_orig.data.num, ROUNDING);
1100 mpfr_set(b, (mpfr_ptr )b_orig.data.num, ROUNDING);
1101 mpfr_set(c, (mpfr_ptr )c_orig.data.num, ROUNDING);
1102 mpfr_set(d, (mpfr_ptr )d_orig.data.num, ROUNDING);
1103 @<Reduce to the case that |a,c>=0|, |b,d>0|@>;
1104 while (1) {
1105 mpfr_div(q,a,d, ROUNDING);
1106 mpfr_div(r,c,b, ROUNDING);
1107 cmp = mpfr_cmp(q,r);
1108 if (cmp) {
1109 if (cmp>1) {
1110 mpfr_set(ret->data.num, one, ROUNDING);
1111 } else {
1112 mpfr_set(ret->data.num, minusone, ROUNDING);
1113 }
1114 goto RETURN;
1115 }
1116 mpfr_remainder(q,a,d, ROUNDING);
1117 mpfr_remainder(r,c,b, ROUNDING);
1118 if (mpfr_zero_p(r)) {
1119 if (mpfr_zero_p(q)) {
1120 mpfr_set(ret->data.num, zero, ROUNDING);
1121 } else {
1122 mpfr_set(ret->data.num, one, ROUNDING);
1123 }
1124 goto RETURN;
1125 }
1126 if (mpfr_zero_p(q)) {
1127 mpfr_set(ret->data.num, minusone, ROUNDING);
1128 goto RETURN;
1129 }
1130 mpfr_set(a,b, ROUNDING);
1131 mpfr_set(b,q, ROUNDING);
1132 mpfr_set(c,d, ROUNDING);
1133 mpfr_set(d,r, ROUNDING);
1134 } /* now |a>d>0| and |c>b>0| */
1135 RETURN:
1136 #if DEBUG
1137 fprintf(stdout, "\n%f = ab_vs_cd(%f,%f,%f,%f)", mp_number_to_double(*ret),
1138 mp_number_to_double(a_orig),mp_number_to_double(b_orig),
1139 mp_number_to_double(c_orig),mp_number_to_double(d_orig));
1140 #endif
1141 mp_check_mpfr_t(mp, ret->data.num);
1142 mpfr_clears(q,r,test,a,b,c,d,(mpfr_ptr)0);
1143 return;
1144 }
1145
1146
1147 @ @<Reduce to the case that |a...@>=
1148 if (mpfr_negative_p(a)) {
1149 mpfr_neg(a, a, ROUNDING);
1150 mpfr_neg(b, b, ROUNDING);
1151 }
1152 if (mpfr_negative_p(c)) {
1153 mpfr_neg(c, c, ROUNDING);
1154 mpfr_neg(d, d, ROUNDING);
1155 }
1156 if (!mpfr_positive_p(d)) {
1157 if (!mpfr_negative_p(b)) {
1158 if ((mpfr_zero_p(a) || mpfr_zero_p(b)) && (mpfr_zero_p(c) || mpfr_zero_p(d)))
1159 mpfr_set(ret->data.num, zero, ROUNDING);
1160 else
1161 mpfr_set(ret->data.num, one, ROUNDING);
1162 goto RETURN;
1163 }
1164 if (mpfr_zero_p(d)) {
1165 if (mpfr_zero_p(a))
1166 mpfr_set(ret->data.num, zero, ROUNDING);
1167 else
1168 mpfr_set(ret->data.num, minusone, ROUNDING);
1169 goto RETURN;
1170 }
1171 mpfr_set(q, a, ROUNDING);
1172 mpfr_set(a, c, ROUNDING);
1173 mpfr_set(c, q, ROUNDING);
1174 mpfr_neg(q, b, ROUNDING);
1175 mpfr_neg(b, d, ROUNDING);
1176 mpfr_set(d, q, ROUNDING);
1177 } else if (!mpfr_positive_p(b)) {
1178 if (mpfr_negative_p(b) && mpfr_positive_p(a)) {
1179 mpfr_set(ret->data.num, minusone, ROUNDING);
1180 goto RETURN;
1181 }
1182 if (mpfr_zero_p(c))
1183 mpfr_set(ret->data.num, zero, ROUNDING);
1184 else
1185 mpfr_set(ret->data.num, minusone, ROUNDING);
1186 goto RETURN;
1187 }
1188
1189 @ Now here's a subroutine that's handy for all sorts of path computations:
1190 Given a quadratic polynomial $B(a,b,c;t)$, the |crossing_point| function
1191 returns the unique |fraction| value |t| between 0 and~1 at which
1192 $B(a,b,c;t)$ changes from positive to negative, or returns
1193 |t=fraction_one+1| if no such value exists. If |a<0| (so that $B(a,b,c;t)$
1194 is already negative at |t=0|), |crossing_point| returns the value zero.
1195
1196 The general bisection method is quite simple when $n=2$, hence
1197 |crossing_point| does not take much time. At each stage in the
1198 recursion we have a subinterval defined by |l| and~|j| such that
1199 $B(a,b,c;2^{-l}(j+t))=B(x_0,x_1,x_2;t)$, and we want to ``zero in'' on
1200 the subinterval where $x_0\G0$ and $\min(x_1,x_2)<0$.
1201
1202 It is convenient for purposes of calculation to combine the values
1203 of |l| and~|j| in a single variable $d=2^l+j$, because the operation
1204 of bisection then corresponds simply to doubling $d$ and possibly
1205 adding~1. Furthermore it proves to be convenient to modify
1206 our previous conventions for bisection slightly, maintaining the
1207 variables $X_0=2^lx_0$, $X_1=2^l(x_0-x_1)$, and $X_2=2^l(x_1-x_2)$.
1208 With these variables the conditions $x_0\ge0$ and $\min(x_1,x_2)<0$ are
1209 equivalent to $\max(X_1,X_1+X_2)>X_0\ge0$.
1210
1211 The following code maintains the invariant relations
1212 $0\L|x0|<\max(|x1|,|x1|+|x2|)$,
1213 $\vert|x1|\vert<2^{30}$, $\vert|x2|\vert<2^{30}$;
1214 it has been constructed in such a way that no arithmetic overflow
1215 will occur if the inputs satisfy
1216 $a<2^{30}$, $\vert a-b\vert<2^{30}$, and $\vert b-c\vert<2^{30}$.
1217
1218 @d no_crossing { mpfr_set(ret->data.num, fraction_one_plus_mpfr_t, ROUNDING); goto RETURN; }
1219 @d one_crossing { mpfr_set(ret->data.num, fraction_one_mpfr_t, ROUNDING); goto RETURN; }
1220 @d zero_crossing { mpfr_set(ret->data.num, zero, ROUNDING); goto RETURN; }
1221
1222 @c
1223 static void mp_binary_crossing_point (MP mp, mp_number *ret, mp_number aa, mp_number bb, mp_number cc) {
1224 mpfr_t a,b,c;
1225 double d; /* recursive counter */
1226 mpfr_t x, xx, x0, x1, x2; /* temporary registers for bisection */
1227 mpfr_t scratch;
1228 mpfr_inits2 (precision_bits, a,b,c, x,xx,x0,x1,x2, scratch,(mpfr_ptr)0);
1229 mpfr_set(a, (mpfr_ptr )aa.data.num, ROUNDING);
1230 mpfr_set(b, (mpfr_ptr )bb.data.num, ROUNDING);
1231 mpfr_set(c, (mpfr_ptr )cc.data.num, ROUNDING);
1232 if (mpfr_negative_p(a))
1233 zero_crossing;
1234 if (!mpfr_negative_p(c)) {
1235 if (!mpfr_negative_p(b)) {
1236 if (mpfr_positive_p(c)) {
1237 no_crossing;
1238 } else if (mpfr_zero_p(a) && mpfr_zero_p(b)) {
1239 no_crossing;
1240 } else {
1241 one_crossing;
1242 }
1243 }
1244 if (mpfr_zero_p(a))
1245 zero_crossing;
1246 } else if (mpfr_zero_p(a)) {
1247 if (!mpfr_positive_p(b))
1248 zero_crossing;
1249 }
1250
1251 /* Use bisection to find the crossing point... */
1252 d = epsilonf;
1253 mpfr_set(x0, a, ROUNDING);
1254 mpfr_sub(x1,a, b, ROUNDING);
1255 mpfr_sub(x2,b, c, ROUNDING);
1256 do {
1257 /* not sure why the error correction has to be >= 1E-12 */
1258 mpfr_add(x, x1, x2, ROUNDING);
1259 mpfr_div(x, x, two_mpfr_t, ROUNDING);
1260 mpfr_add_d (x, x, 1E-12, ROUNDING);
1261 mpfr_sub(scratch, x1, x0, ROUNDING);
1262 if (mpfr_greater_p(scratch, x0)) {
1263 mpfr_set(x2, x, ROUNDING);
1264 mpfr_add(x0, x0, x0, ROUNDING);
1265 d += d;
1266 } else {
1267 mpfr_add(xx, scratch, x, ROUNDING);
1268 if (mpfr_greater_p(xx,x0)) {
1269 mpfr_set(x2,x, ROUNDING);
1270 mpfr_add(x0, x0, x0, ROUNDING);
1271 d += d;
1272 } else {
1273 mpfr_sub(x0, x0, xx, ROUNDING);
1274 if (!mpfr_greater_p(x,x0)) {
1275 mpfr_add(scratch, x, x2, ROUNDING);
1276 if (!mpfr_greater_p(scratch, x0))
1277 no_crossing;
1278 }
1279 mpfr_set(x1,x, ROUNDING);
1280 d = d + d + epsilonf;
1281 }
1282 }
1283 } while (d < fraction_one);
1284 mpfr_set_d(scratch, d, ROUNDING);
1285 mpfr_sub(ret->data.num,scratch, fraction_one_mpfr_t, ROUNDING);
1286 RETURN:
1287 #if DEBUG
1288 fprintf(stdout, "\n%f = crossing_point(%f,%f,%f)", mp_number_to_double(*ret),
1289 mp_number_to_double(aa),mp_number_to_double(bb),mp_number_to_double(cc));
1290 #endif
1291 mpfr_clears (a,b,c, x,xx,x0,x1,x2, scratch, (mpfr_ptr)0);
1292 mp_check_mpfr_t(mp, ret->data.num);
1293 return;
1294 }
1295
1296
1297 @ We conclude this set of elementary routines with some simple rounding
1298 and truncation operations.
1299
1300
1301 @ |round_unscaled| rounds a |scaled| and converts it to |int|
1302 @c
1303 int mp_round_unscaled(mp_number x_orig) {
1304 double xx = mp_number_to_double(x_orig);
1305 int x = (int)ROUND(xx);
1306 return x;
1307 }
1308
1309 @ |number_floor| floors a number
1310
1311 @c
1312 void mp_number_floor (mp_number *i) {
1313 mpfr_rint_floor(i->data.num, i->data.num, MPFR_RNDD);
1314 }
1315
1316 @ |fraction_to_scaled| rounds a |fraction| and converts it to |scaled|
1317 @c
1318 void mp_binary_fraction_to_round_scaled (mp_number *x_orig) {
1319 x_orig->type = mp_scaled_type;
1320 mpfr_div(x_orig->data.num, x_orig->data.num, fraction_multiplier_mpfr_t, ROUNDING);
1321 }
1322
1323
1324
1325 @* Algebraic and transcendental functions.
1326 \MP\ computes all of the necessary special functions from scratch, without
1327 relying on |real| arithmetic or system subroutines for sines, cosines, etc.
1328
1329 @
1330
1331 @c
1332 void mp_binary_square_rt (MP mp, mp_number *ret, mp_number x_orig) { /* return, x: scaled */
1333 if (!mpfr_positive_p((mpfr_ptr)x_orig.data.num)) {
1334 @<Handle square root of zero or negative argument@>;
1335 } else {
1336 mpfr_sqrt(ret->data.num, x_orig.data.num, ROUNDING);
1337 }
1338 mp_check_mpfr_t(mp, ret->data.num);
1339 }
1340
1341
1342 @ @<Handle square root of zero...@>=
1343 {
1344 if (mpfr_negative_p((mpfr_ptr)x_orig.data.num)) {
1345 char msg[256];
1346 const char *hlp[] = {
1347 "Since I don't take square roots of negative numbers,",
1348 "I'm zeroing this one. Proceed, with fingers crossed.",
1349 NULL };
1350 char *xstr = mp_binary_number_tostring (mp, x_orig);
1351 mp_snprintf(msg, 256, "Square root of %s has been replaced by 0", xstr);
1352 free(xstr);
1353 @.Square root...replaced by 0@>;
1354 mp_error (mp, msg, hlp, true);
1355 }
1356 mpfr_set_zero(ret->data.num,1); /* 1 == positive */
1357 return;
1358 }
1359
1360
1361 @ Pythagorean addition $\psqrt{a^2+b^2}$ is implemented by a quick hack
1362
1363 @c
1364 void mp_binary_pyth_add (MP mp, mp_number *ret, mp_number a_orig, mp_number b_orig) {
1365 mpfr_t a, b, asq, bsq;
1366 mpfr_inits2(precision_bits, a,b, asq, bsq, (mpfr_ptr)0);
1367 mpfr_set(a, (mpfr_ptr)a_orig.data.num, ROUNDING);
1368 mpfr_set(b, (mpfr_ptr)b_orig.data.num, ROUNDING);
1369 mpfr_mul(asq, a, a, ROUNDING);
1370 mpfr_mul(bsq, b, b, ROUNDING);
1371 mpfr_add(a, asq, bsq, ROUNDING);
1372 mpfr_sqrt(ret->data.num, a, ROUNDING);
1373 mp_check_mpfr_t(mp, ret->data.num);
1374 mpfr_clears(a,b, asq, bsq, (mpfr_ptr)0);
1375 }
1376
1377 @ Here is a similar algorithm for $\psqrt{a^2-b^2}$. Same quick hack, also.
1378
1379 @c
1380 void mp_binary_pyth_sub (MP mp, mp_number *ret, mp_number a_orig, mp_number b_orig) {
1381 mpfr_t a, b, asq, bsq;
1382 mpfr_inits2(precision_bits, a,b, asq, bsq, (mpfr_ptr)0);
1383 mpfr_set(a, (mpfr_ptr)a_orig.data.num, ROUNDING);
1384 mpfr_set(b, (mpfr_ptr)b_orig.data.num, ROUNDING);
1385 if (!mpfr_greater_p(a,b)) {
1386 @<Handle erroneous |pyth_sub| and set |a:=0|@>;
1387 } else {
1388 mpfr_mul(asq, a, a, ROUNDING);
1389 mpfr_mul(bsq, b, b, ROUNDING);
1390 mpfr_sub(a, asq, bsq, ROUNDING);
1391 mpfr_sqrt(a, a, ROUNDING);
1392 }
1393 mpfr_set(ret->data.num, a, ROUNDING);
1394 mp_check_mpfr_t(mp, ret->data.num);
1395 }
1396
1397
1398 @ @<Handle erroneous |pyth_sub| and set |a:=0|@>=
1399 {
1400 if (mpfr_less_p(a, b)) {
1401 char msg[256];
1402 const char *hlp[] = {
1403 "Since I don't take square roots of negative numbers,",
1404 "I'm zeroing this one. Proceed, with fingers crossed.",
1405 NULL };
1406 char *astr = mp_binary_number_tostring (mp, a_orig);
1407 char *bstr = mp_binary_number_tostring (mp, b_orig);
1408 mp_snprintf (msg, 256, "Pythagorean subtraction %s+-+%s has been replaced by 0", astr, bstr);
1409 free(astr);
1410 free(bstr);
1411 @.Pythagorean...@>;
1412 mp_error (mp, msg, hlp, true);
1413 }
1414 mpfr_set_zero(a,1); /* 1 == positive */
1415 }
1416
1417
1418 @ Here is the routine that calculates $2^8$ times the natural logarithm
1419 of a |scaled| quantity;
1420
1421 @c
1422 void mp_binary_m_log (MP mp, mp_number *ret, mp_number x_orig) {
1423 if (!mpfr_positive_p((mpfr_ptr)x_orig.data.num)) {
1424 @<Handle non-positive logarithm@>;
1425 } else {
1426 mpfr_log(ret->data.num, x_orig.data.num, ROUNDING);
1427 mp_check_mpfr_t(mp, ret->data.num);
1428 mpfr_mul_si(ret->data.num, ret->data.num, 256, ROUNDING);
1429 }
1430 mp_check_mpfr_t(mp, ret->data.num);
1431 }
1432
1433 @ @<Handle non-positive logarithm@>=
1434 {
1435 char msg[256];
1436 const char *hlp[] = {
1437 "Since I don't take logs of non-positive numbers,",
1438 "I'm zeroing this one. Proceed, with fingers crossed.",
1439 NULL };
1440 char *xstr = mp_binary_number_tostring (mp, x_orig);
1441 mp_snprintf (msg, 256, "Logarithm of %s has been replaced by 0", xstr);
1442 free (xstr);
1443 @.Logarithm...replaced by 0@>;
1444 mp_error (mp, msg, hlp, true);
1445 mpfr_set_zero(ret->data.num,1); /* 1 == positive */
1446 }
1447
1448
1449 @ Conversely, the exponential routine calculates $\exp(x/2^8)$,
1450 when |x| is |scaled|.
1451
1452 @c
1453 void mp_binary_m_exp (MP mp, mp_number *ret, mp_number x_orig) {
1454 mpfr_t temp;
1455 mpfr_init2(temp, precision_bits);
1456 mpfr_div_si(temp, x_orig.data.num, 256, ROUNDING);
1457 mpfr_exp(ret->data.num, temp, ROUNDING);
1458 mp_check_mpfr_t(mp, ret->data.num);
1459 mpfr_clear (temp);
1460 }
1461
1462
1463 @ Given integers |x| and |y|, not both zero, the |n_arg| function
1464 returns the |angle| whose tangent points in the direction $(x,y)$.
1465
1466 @c
1467 void mp_binary_n_arg (MP mp, mp_number *ret, mp_number x_orig, mp_number y_orig) {
1468 if (mpfr_zero_p((mpfr_ptr )x_orig.data.num) && mpfr_zero_p((mpfr_ptr )y_orig.data.num)) {
1469 @<Handle undefined arg@>;
1470 } else {
1471 mpfr_t atan2val, oneeighty_angle;
1472 mpfr_init2(atan2val, precision_bits);
1473 mpfr_init2(oneeighty_angle, precision_bits);
1474 ret->type = mp_angle_type;
1475 mpfr_set_si(oneeighty_angle, 180 * angle_multiplier, ROUNDING);
1476 mpfr_div(oneeighty_angle, oneeighty_angle, PI_mpfr_t, ROUNDING);
1477 checkZero((mpfr_ptr)y_orig.data.num);
1478 checkZero((mpfr_ptr)x_orig.data.num);
1479 mpfr_atan2(atan2val, y_orig.data.num, x_orig.data.num, ROUNDING);
1480 mpfr_mul(ret->data.num, atan2val, oneeighty_angle, ROUNDING);
1481 checkZero((mpfr_ptr)ret->data.num);
1482 mpfr_clear(atan2val);
1483 mpfr_clear(oneeighty_angle);
1484 }
1485 mp_check_mpfr_t(mp, ret->data.num);
1486 }
1487
1488
1489 @ @<Handle undefined arg@>=
1490 {
1491 const char *hlp[] = {
1492 "The `angle' between two identical points is undefined.",
1493 "I'm zeroing this one. Proceed, with fingers crossed.",
1494 NULL };
1495 mp_error (mp, "angle(0,0) is taken as zero", hlp, true);
1496 @.angle(0,0)...zero@>;
1497 mpfr_set_zero(ret->data.num,1); /* 1 == positive */
1498 }
1499
1500
1501 @ Conversely, the |n_sin_cos| routine takes an |angle| and produces the sine
1502 and cosine of that angle. The results of this routine are
1503 stored in global integer variables |n_sin| and |n_cos|.
1504
1505 @ Calculate sines and cosines.
1506
1507 @c
1508 void mp_binary_sin_cos (MP mp, mp_number z_orig, mp_number *n_cos, mp_number *n_sin) {
1509 mpfr_t rad;
1510 mpfr_t one_eighty;
1511 mpfr_init2(rad, precision_bits);
1512 mpfr_init2(one_eighty, precision_bits);
1513 mpfr_set_si(one_eighty, 180 * 16, ROUNDING);
1514 mpfr_mul (rad, z_orig.data.num, PI_mpfr_t, ROUNDING);
1515 mpfr_div (rad, rad, one_eighty, ROUNDING);
1516
1517 mpfr_sin (n_sin->data.num, rad, ROUNDING);
1518 mpfr_cos (n_cos->data.num, rad, ROUNDING);
1519
1520 mpfr_mul (n_cos->data.num,n_cos->data.num, fraction_multiplier_mpfr_t, ROUNDING);
1521 mpfr_mul (n_sin->data.num,n_sin->data.num, fraction_multiplier_mpfr_t, ROUNDING);
1522 mp_check_mpfr_t(mp, n_cos->data.num);
1523 mp_check_mpfr_t(mp, n_sin->data.num);
1524 mpfr_clear (rad);
1525 mpfr_clear (one_eighty);
1526 }
1527
1528 @ This is the http://www-cs-faculty.stanford.edu/~uno/programs/rng.c
1529 with small cosmetic modifications.
1530
1531 @c
1532 #define KK 100 /* the long lag */
1533 #define LL 37 /* the short lag */
1534 #define MM (1L<<30) /* the modulus */
1535 #define mod_diff(x,y) (((x)-(y))&(MM-1)) /* subtraction mod MM */
1536 /* */
1537 static long ran_x[KK]; /* the generator state */
1538 /* */
1539 static void ran_array(long aa[],int n) /* put n new random numbers in aa */
1540 /* long aa[] destination */
1541 /* int n array length (must be at least KK) */
1542 {
1543 register int i,j;
1544 for (j=0;j<KK;j++) aa[j]=ran_x[j];
1545 for (;j<n;j++) aa[j]=mod_diff(aa[j-KK],aa[j-LL]);
1546 for (i=0;i<LL;i++,j++) ran_x[i]=mod_diff(aa[j-KK],aa[j-LL]);
1547 for (;i<KK;i++,j++) ran_x[i]=mod_diff(aa[j-KK],ran_x[i-LL]);
1548 }
1549 /* */
1550 /* the following routines are from exercise 3.6--15 */
1551 /* after calling ran_start, get new randoms by, e.g., "x=ran_arr_next()" */
1552 /* */
1553 #define QUALITY 1009 /* recommended quality level for high-res use */
1554 static long ran_arr_buf[QUALITY];
1555 static long ran_arr_dummy=-1, ran_arr_started=-1;
1556 static long *ran_arr_ptr=&ran_arr_dummy; /* the next random number, or -1 */
1557 /* */
1558 #define TT 70 /* guaranteed separation between streams */
1559 #define is_odd(x) ((x)&1) /* units bit of x */
1560 /* */
1561 static void ran_start(long seed) /* do this before using ran_array */
1562 /* long seed selector for different streams */
1563 {
1564 register int t,j;
1565 long x[KK+KK-1]; /* the preparation buffer */
1566 register long ss=(seed+2)&(MM-2);
1567 for (j=0;j<KK;j++) {
1568 x[j]=ss; /* bootstrap the buffer */
1569 ss<<=1; if (ss>=MM) ss-=MM-2; /* cyclic shift 29 bits */
1570 }
1571 x[1]++; /* make x[1] (and only x[1]) odd */
1572 for (ss=seed&(MM-1),t=TT-1; t; ) {
1573 for (j=KK-1;j>0;j--) x[j+j]=x[j], x[j+j-1]=0; /* "square" */
1574 for (j=KK+KK-2;j>=KK;j--)
1575 x[j-(KK-LL)]=mod_diff(x[j-(KK-LL)],x[j]),
1576 x[j-KK]=mod_diff(x[j-KK],x[j]);
1577 if (is_odd(ss)) { /* "multiply by z" */
1578 for (j=KK;j>0;j--) x[j]=x[j-1];
1579 x[0]=x[KK]; /* shift the buffer cyclically */
1580 x[LL]=mod_diff(x[LL],x[KK]);
1581 }
1582 if (ss) ss>>=1; else t--;
1583 }
1584 for (j=0;j<LL;j++) ran_x[j+KK-LL]=x[j];
1585 for (;j<KK;j++) ran_x[j-LL]=x[j];
1586 for (j=0;j<10;j++) ran_array(x,KK+KK-1); /* warm things up */
1587 ran_arr_ptr=&ran_arr_started;
1588 }
1589 /* */
1590 #define ran_arr_next() (*ran_arr_ptr>=0? *ran_arr_ptr++: ran_arr_cycle())
1591 static long ran_arr_cycle(void)
1592 {
1593 if (ran_arr_ptr==&ran_arr_dummy)
1594 ran_start(314159L); /* the user forgot to initialize */
1595 ran_array(ran_arr_buf,QUALITY);
1596 ran_arr_buf[KK]=-1;
1597 ran_arr_ptr=ran_arr_buf+1;
1598 return ran_arr_buf[0];
1599 }
1600
1601
1602
1603
1604 @ To initialize the |randoms| table, we call the following routine.
1605
1606 @c
1607 void mp_init_randoms (MP mp, int seed) {
1608 int j, jj, k; /* more or less random integers */
1609 int i; /* index into |randoms| */
1610 j = abs (seed);
1611 while (j >= fraction_one) {
1612 j = j/2;
1613 }
1614 k = 1;
1615 for (i = 0; i <= 54; i++) {
1616 jj = k;
1617 k = j - k;
1618 j = jj;
1619 if (k<0)
1620 k += fraction_one;
1621 mpfr_set_si(mp->randoms[(i * 21) % 55].data.num, j, ROUNDING);
1622 }
1623 mp_new_randoms (mp);
1624 mp_new_randoms (mp);
1625 mp_new_randoms (mp); /* ``warm up'' the array */
1626
1627 ran_start ((unsigned long)seed);
1628
1629 }
1630
1631 @ @c
1632 void mp_binary_number_modulo (mp_number *a, mp_number b) {
1633 mpfr_remainder (a->data.num, a->data.num, b.data.num, ROUNDING);
1634 }
1635
1636 @ To consume a random integer for the uniform generator, the program below will say `|next_unif_random|'.
1637
1638 @c
1639 static void mp_next_unif_random (MP mp, mp_number *ret) {
1640 mp_number rop;
1641 unsigned long int op;
1642 float flt_op ;
1643 (void)mp;
1644 mp_new_number (mp, &rop, mp_scaled_type);
1645 op = (unsigned)ran_arr_next();
1646 flt_op = op/(MM*1.0);
1647 mpfr_set_d ((mpfr_ptr)(rop.data.num), flt_op,ROUNDING);
1648 mp_number_clone (ret, rop);
1649 free_number (rop);
1650 }
1651
1652
1653
1654 @ To consume a random fraction, the program below will say `|next_random|'.
1655
1656 @c
1657 static void mp_next_random (MP mp, mp_number *ret) {
1658 if ( mp->j_random==0 )
1659 mp_new_randoms(mp);
1660 else
1661 mp->j_random = mp->j_random-1;
1662 mp_number_clone (ret, mp->randoms[mp->j_random]);
1663 }
1664
1665 @ To produce a uniform random number in the range |0<=u<x| or |0>=u>x|
1666 or |0=u=x|, given a |scaled| value~|x|, we proceed as shown here.
1667
1668 Note that the call of |take_fraction| will produce the values 0 and~|x|
1669 with about half the probability that it will produce any other particular
1670 values between 0 and~|x|, because it rounds its answers.
1671
1672 @c
1673 static void mp_binary_m_unif_rand (MP mp, mp_number *ret, mp_number x_orig) {
1674 mp_number y; /* trial value */
1675 mp_number x, abs_x;
1676 mp_number u;
1677 char *r ;mpfr_exp_t e;
1678 new_fraction (y);
1679 new_number (x);
1680 new_number (abs_x);
1681 new_number (u);
1682 mp_number_clone (&x, x_orig);
1683 mp_number_clone (&abs_x, x);
1684 mp_binary_abs (&abs_x);
1685 mp_next_unif_random(mp, &u);
1686 mpfr_mul (y.data.num, abs_x.data.num, u.data.num, ROUNDING);
1687 free_number (u);
1688 if (mp_number_equal(y, abs_x)) {
1689 mp_number_clone (ret, ((math_data *)mp->math)->zero_t);
1690 } else if (mp_number_greater(x, ((math_data *)mp->math)->zero_t)) {
1691 mp_number_clone (ret, y);
1692 } else {
1693 mp_number_clone (ret, y);
1694 mp_number_negate (ret);
1695 }
1696 r = mpfr_get_str(NULL, // char *str,
1697 &e, // mpfr_exp_t *expptr,
1698 10, // int b,
1699 0, // size_t n,
1700 ret->data.num, // mpfr_t op,
1701 ROUNDING // mpfr_rnd_t rnd
1702 );
1703 printf("\nret=%s e=%ld\n",r,e);
1704 mpfr_free_str(r);
1705 free_number (abs_x);
1706 free_number (x);
1707 free_number (y);
1708 }
1709
1710
1711
1712 @ Finally, a normal deviate with mean zero and unit standard deviation
1713 can readily be obtained with the ratio method (Algorithm 3.4.1R in
1714 {\sl The Art of Computer Programming\/}).
1715
1716 @c
1717 static void mp_binary_m_norm_rand (MP mp, mp_number *ret) {
1718 mp_number ab_vs_cd;
1719 mp_number abs_x;
1720 mp_number u;
1721 mp_number r;
1722 mp_number la, xa;
1723 new_number (ab_vs_cd);
1724 new_number (la);
1725 new_number (xa);
1726 new_number (abs_x);
1727 new_number (u);
1728 new_number (r);
1729
1730 do {
1731 do {
1732 mp_number v;
1733 new_number (v);
1734 mp_next_random(mp, &v);
1735 mp_number_substract (&v, ((math_data *)mp->math)->fraction_half_t);
1736 mp_binary_number_take_fraction (mp,&xa, ((math_data *)mp->math)->sqrt_8_e_k, v);
1737 free_number (v);
1738 mp_next_random(mp, &u);
1739 mp_number_clone (&abs_x, xa);
1740 mp_binary_abs (&abs_x);
1741 } while (!mp_number_less(abs_x, u));
1742 mp_binary_number_make_fraction (mp, &r, xa, u);
1743 mp_number_clone (&xa, r);
1744 mp_binary_m_log (mp,&la, u);
1745 mp_set_binary_from_substraction(&la, ((math_data *)mp->math)->twelve_ln_2_k, la);
1746 mp_binary_ab_vs_cd (mp,&ab_vs_cd, ((math_data *)mp->math)->one_k, la, xa, xa);
1747 /*mp_ab_vs_cd (mp,&ab_vs_cd, ((math_data *)mp->math)->one_k, la, xa, xa);*/
1748 } while (mp_number_less(ab_vs_cd,((math_data *)mp->math)->zero_t));
1749 mp_number_clone (ret, xa);
1750 free_number (ab_vs_cd);
1751 free_number (r);
1752 free_number (abs_x);
1753 free_number (la);
1754 free_number (xa);
1755 free_number (u);
1756 }
1757
1758
1759
1760 @ The following subroutine is used only in norm_rand and tests if $ab$ is
1761 greater than, equal to, or less than~$cd$.
1762 The result is $+1$, 0, or~$-1$ in the three respective cases.
1763
1764 @c
1765 static void mp_binary_ab_vs_cd (MP mp, mp_number *ret, mp_number a_orig, mp_number b_orig, mp_number c_orig, mp_number d_orig) {
1766 mpfr_t a, b, c, d;
1767 mpfr_t ab, cd;
1768
1769 int cmp = 0;
1770 (void)mp;
1771 mpfr_inits2(precision_bits, a,b,c,d,ab,cd,(mpfr_ptr)0);
1772 mpfr_set(a, (mpfr_ptr )a_orig.data.num, ROUNDING);
1773 mpfr_set(b, (mpfr_ptr )b_orig.data.num, ROUNDING);
1774 mpfr_set(c, (mpfr_ptr )c_orig.data.num, ROUNDING);
1775 mpfr_set(d, (mpfr_ptr )d_orig.data.num, ROUNDING);
1776
1777 mpfr_mul(ab,a,b, ROUNDING);
1778 mpfr_mul(cd,c,d, ROUNDING);
1779
1780 mpfr_set(ret->data.num, zero, ROUNDING);
1781 cmp = mpfr_cmp(ab,cd);
1782 if (cmp) {
1783 if (cmp>0)
1784 mpfr_set(ret->data.num, one, ROUNDING);
1785 else
1786 mpfr_set(ret->data.num, minusone, ROUNDING);
1787 }
1788 mp_check_mpfr_t(mp, ret->data.num);
1789 mpfr_clears(a,b,c,d,ab,cd,(mpfr_ptr)0);
1790 return;
1791 }
1792
1793
LuaTeX git mirror