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