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.
7 \def\title
{Math support functions for MPFR based math
}
13 #include
<w2c
/config.h
>
18 #include
"mpmathbinary.h" /* internal header
*/
19 #define
ROUND(a
) floor
((a
)+0.5)
26 #ifndef MPMATHBINARY_H
27 #define MPMATHBINARY_H
1
29 #include
"mpmp.h" /* internal header
*/
32 @
<Internal library declarations@
>;
35 @
* Math initialization.
37 First
, here are some very important constants.
40 @d E_STRING
"2.7182818284590452353602874713526624977572470936999595749669676277240766303535"
41 @d PI_STRING
"3.1415926535897932384626433832795028841971693993751058209749445923078164062862"
42 @d fraction_multiplier
4096
43 @d angle_multiplier
16
45 @ Here are the functions that are static as they are not used elsewhere
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
);
124 @ We do not want special numbers as return values for functions
, so
:
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);
133 int binary_number_check
(mpfr_t dec
)
136 if
(!mpfr_number_p
(dec
)) {
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
);
144 mpfr_set_zero
(dec
,1); /* 1 == positive
*/
150 void mp_check_mpfr_t
(MP mp
, mpfr_t dec
)
152 mp-
>arith_error
= binary_number_check
(dec
);
158 @ Precision IO uses |double| because |MPFR_PREC_MAX| overflows int.
161 static double precision_bits
;
162 mpfr_prec_t precision_digits_to_bits
(double i
)
166 double precision_bits_to_digits
(mpfr_prec_t d
)
172 @ And these are the ones that
{\it are
} used elsewhere
174 @
<Internal library declarations@
>=
175 void
* mp_initialize_binary_math
(MP mp
);
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
*/
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)
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
;
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
);
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);
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.
245 @d MAX_PRECISION
1000.0
246 @d DEF_PRECISION
34.0
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
();
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);
285 mp_new_number
(mp
, &math->arc_tol_k, mp_fraction_type);
287 mpfr_div_si
(math-
>arc_tol_k.data.num
, one
, 4096, ROUNDING
);
288 /* quit when change in arc length estimate reaches this
*/
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
);
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);
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$
*/
311 mp_new_number
(mp
, &math->twelve_ln_2_k, mp_fraction_type);
313 mpfr_set_d
(math-
>twelve_ln_2_k.data.num
, 139548959.6165 / 65536.0, ROUNDING
);
314 /* $
2^
{24}\cdot12\ln2\approx139548959.6165$
*/
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);
322 mpfr_set_d
(math-
>twelvebits_3.data.num
, 1365 / 65536.0, ROUNDING
);
323 /* $
1365\approx
2^
{12}/3$
*/
325 mp_new_number
(mp
, &math->twentysixbits_sqrt2_t, mp_fraction_type);
327 mpfr_set_d
(math-
>twentysixbits_sqrt2_t.data.num
, 94906265.62 / 65536.0, ROUNDING
);
328 /* $
2^
{26}\sqrt2\approx94906265.62$
*/
330 mp_new_number
(mp
, &math->twentyeightbits_d_t, mp_fraction_type);
332 mpfr_set_d
(math-
>twentyeightbits_d_t.data.num
, 35596754.69 / 65536.0, ROUNDING
);
333 /* $
2^
{28}d\approx35596754.69$
*/
335 mp_new_number
(mp
, &math->twentysevenbits_sqrt2_d_t, mp_fraction_type);
337 mpfr_set_d
(math-
>twentysevenbits_sqrt2_d_t.data.num
, 25170706.63 / 65536.0, ROUNDING
);
338 /* $
2^
{27}\sqrt2\
,d\approx25170706.63$
*/
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
);
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
;
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
);
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
();
460 @ Creating an destroying |mp_number| objects
463 void mp_new_number
(MP mp
, mp_number
*n
, mp_number_type t
) {
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
*/
474 void mp_free_number
(MP mp
, mp_number
*n
) {
477 mpfr_clear
(n-
>data.num
);
480 n-
>type
= mp_nan_type
;
483 @ Here are the low-level functions on |mp_number| items
, setters first.
486 void mp_set_binary_from_int
(mp_number
*A
, int B
) {
487 mpfr_set_si
(A-
>data.num
,B
, ROUNDING
);
489 void mp_set_binary_from_boolean
(mp_number
*A
, int B
) {
490 mpfr_set_si
(A-
>data.num
,B
, ROUNDING
);
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
);
496 void mp_set_binary_from_double
(mp_number
*A
, double B
) {
497 mpfr_set_d
(A-
>data.num
, B
, ROUNDING
);
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
);
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
);
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
);
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
);
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
);
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
);
517 void mp_set_binary_from_of_the_way
(MP mp
, mp_number
*A
, mp_number t
, mp_number B
, mp_number C
) {
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
);
526 mp_check_mpfr_t
(mp
, A-
>data.num
);
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
);
532 void mp_number_add
(mp_number
*A
, mp_number B
) {
533 mpfr_add
(A-
>data.num
,A-
>data.num
,B.data.num
, ROUNDING
);
535 void mp_number_substract
(mp_number
*A
, mp_number B
) {
536 mpfr_sub
(A-
>data.num
,A-
>data.num
,B.data.num
, ROUNDING
);
538 void mp_number_half
(mp_number
*A
) {
539 mpfr_div_si
(A-
>data.num
, A-
>data.num
, 2, ROUNDING
);
541 void mp_number_halfp
(mp_number
*A
) {
542 mpfr_div_si
(A-
>data.num
,A-
>data.num
, 2, ROUNDING
);
544 void mp_number_double
(mp_number
*A
) {
545 mpfr_mul_si
(A-
>data.num
,A-
>data.num
, 2, ROUNDING
);
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
);
550 void mp_number_multiply_int
(mp_number
*A
, int B
) {
551 mpfr_mul_si
(A-
>data.num
,A-
>data.num
, B
, ROUNDING
);
553 void mp_number_divide_int
(mp_number
*A
, int B
) {
554 mpfr_div_si
(A-
>data.num
,A-
>data.num
, B
, ROUNDING
);
556 void mp_binary_abs
(mp_number
*A
) {
557 mpfr_abs
(A-
>data.num
, A-
>data.num
, ROUNDING
);
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
);
563 void mp_number_swap
(mp_number
*A
, mp_number
*B
) {
564 mpfr_swap
(A-
>data.num
, B-
>data.num
);
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
);
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
);
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
);
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
);
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
!
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);
598 @d odd
(A
) (abs
(A
)%2==1)
601 int mp_number_to_int
(mp_number A
) {
603 if
(mpfr_fits_sint_p
(A.data.num
, ROUNDING
)) {
604 result
= mpfr_get_si
(A.data.num
, ROUNDING
);
608 int mp_number_to_boolean
(mp_number A
) {
610 if
(mpfr_fits_sint_p
(A.data.num
, ROUNDING
)) {
611 result
= mpfr_get_si
(A.data.num
, ROUNDING
);
615 double mp_number_to_double
(mp_number A
) {
617 if
(mpfr_number_p
(A.data.num
)) {
618 res
= mpfr_get_d
(A.data.num
, ROUNDING
);
622 int mp_number_odd
(mp_number A
) {
623 return odd
(mp_number_to_int
(A
));
625 int mp_number_equal
(mp_number A
, mp_number B
) {
626 return mpfr_equal_p
(A.data.num
,B.data.num
);
628 int mp_number_greater
(mp_number A
, mp_number B
) {
629 return mpfr_greater_p
(A.data.num
,B.data.num
);
631 int mp_number_less
(mp_number A
, mp_number B
) {
632 return mpfr_less_p
(A.data.num
,B.data.num
);
634 int mp_number_nonequalabs
(mp_number A
, mp_number B
) {
635 return
!(mpfr_cmpabs
(A.data.num
, B.data.num
)==0);
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.
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|'.
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
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.
661 char
* mp_binnumber_tostring
(mpfr_t n
) {
662 char
*str
= NULL, *buffer
= NULL;
665 if
((str
= mpfr_get_str
(NULL, &exp, 10, 0, n, ROUNDING))>0) {
666 int numprecdigits
= precision_bits_to_digits
(precision_bits
);
670 while
(strlen
(str
)>0 && *(str+strlen(str)-1) == '0' ) {
671 *(str
+strlen
(str
)-1) = '\
0'
; /* get rid of trailing zeroes
*/
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
685 if
(strlen
(str
+j
) == 0) {
689 if
(exp
<=numprecdigits
&& exp > -6) {
691 buffer
[i
++] = str
[j
++];
693 buffer
[i
++] = (str
[j
] ? str
[j
++] : '
0'
);
698 buffer
[i
++] = str
[j
++];
706 while
(absexp--
> 0) {
710 buffer
[i
++] = str
[j
++];
714 buffer
[i
++] = str
[j
++];
718 buffer
[i
++] = str
[j
++];
724 mp_snprintf
(msg
, 256, "%s%d", (exp
>0?
"+":""), (int
)(exp
>0 ?
(exp-1
) : (exp-1
)));
727 buffer
[i
++] = msg
[k
++];
738 char
* mp_binary_number_tostring
(MP mp
, mp_number n
) {
739 return mp_binnumber_tostring
(n.data.num
);
744 void mp_binary_print_number
(MP mp
, mp_number n
) {
745 char
*str
= mp_binary_number_tostring
(mp
, n
);
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
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
);
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.
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.
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.
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@
>
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.
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
);
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
);
811 void mp_binary_make_fraction
(MP mp
, mpfr_t ret
, mpfr_t p
, mpfr_t q
);
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
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@
>
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
);
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
);
833 void mp_binary_take_fraction
(MP mp
, mpfr_t ret
, mpfr_t p
, mpfr_t q
);
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$.
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.
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
);
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.
)
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
);
864 @d halfp
(A
) (integer
)((unsigned
)(A
) >> 1)
866 @
* Scanning numbers in the input
868 The definitions below are temporarily here
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
)
874 static void mp_wrapup_numeric_token
(MP mp
, unsigned char
*start
, unsigned char
*stop
);
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
)
885 void mp_wrapup_numeric_token
(MP mp
, unsigned char
*start
, unsigned char
*stop
) {
888 size_t l
= stop-start
+1;
889 unsigned long lp
, lpbit
;
891 char
*buf
= mp_xmalloc
(mp
, l
+1, 1);
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
++;}
902 lp
= strchr
(bufp
,'.'
) ? lp-1
: lp
;
903 /* strip also trailing
0s
*/
905 while
(*bufp
== '
0'
) {bufp--
; lp
=( ((lp
==0)||
(lp
==1))?
1:lp-1
);}
906 /* at least one digit
, even if the number is
0 */
908 /* bits needed for buf
*/
909 lpbit
= (unsigned long
)ceil
(lp
/log10
(2)+1);
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
)) {
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.)",
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
);
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:",
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
));
938 set_cur_cmd
((mp_variable_type
)mp_numeric_token
);
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--
;
953 if
(mp-
>buffer
[mp-
>cur_input.loc_field
] == '
+' ||
954 mp-
>buffer
[mp-
>cur_input.loc_field
] == '
-'
) {
955 mp-
>cur_input.loc_field
++;
957 while
(mp-
>char_class
[mp-
>buffer
[mp-
>cur_input.loc_field
]] == digit_class
) {
958 mp-
>cur_input.loc_field
++;
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];
965 while
(mp-
>char_class
[mp-
>buffer
[mp-
>cur_input.loc_field
]] == digit_class
) {
966 mp-
>cur_input.loc_field
++;
969 stop
= &mp->buffer[mp->cur_input.loc_field-1];
970 mp_wrapup_numeric_token
(mp
, start
, stop
);
974 @ We just have to collect bytes.
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];
980 while
(mp-
>char_class
[mp-
>buffer
[mp-
>cur_input.loc_field
]] == digit_class
) {
981 mp-
>cur_input.loc_field
++;
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
++;
991 stop
= &mp->buffer[mp->cur_input.loc_field-1];
992 mp_wrapup_numeric_token
(mp
, start
, stop
);
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
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
)
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.
)
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.
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.
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
*/
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
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
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
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
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
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
1072 if
(!mpfr_equal_p
(t.data.num
, one
)) { // t
!= 1
1073 mpfr_div
(num
, num
, t.data.num
, ROUNDING
); // num
= num
/ t
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
);
1080 mp_binary_make_fraction
(mp
, ret-
>data.num
, num
, denom
);
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
);
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.
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
*/
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|@
>;
1105 mpfr_div
(q
,a
,d
, ROUNDING
);
1106 mpfr_div
(r
,c
,b
, ROUNDING
);
1107 cmp
= mpfr_cmp
(q
,r
);
1110 mpfr_set
(ret-
>data.num
, one
, ROUNDING
);
1112 mpfr_set
(ret-
>data.num
, minusone
, ROUNDING
);
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
);
1122 mpfr_set
(ret-
>data.num
, one
, ROUNDING
);
1126 if
(mpfr_zero_p
(q
)) {
1127 mpfr_set
(ret-
>data.num
, minusone
, ROUNDING
);
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|
*/
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
));
1141 mp_check_mpfr_t
(mp
, ret-
>data.num
);
1142 mpfr_clears
(q
,r
,test
,a
,b
,c
,d
,(mpfr_ptr
)0);
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
);
1152 if
(mpfr_negative_p
(c
)) {
1153 mpfr_neg
(c
, c
, ROUNDING
);
1154 mpfr_neg
(d
, d
, ROUNDING
);
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
);
1161 mpfr_set
(ret-
>data.num
, one
, ROUNDING
);
1164 if
(mpfr_zero_p
(d
)) {
1166 mpfr_set
(ret-
>data.num
, zero
, ROUNDING
);
1168 mpfr_set
(ret-
>data.num
, minusone
, ROUNDING
);
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
);
1183 mpfr_set
(ret-
>data.num
, zero
, ROUNDING
);
1185 mpfr_set
(ret-
>data.num
, minusone
, ROUNDING
);
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.
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$.
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$.
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}$.
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; }
1223 static void mp_binary_crossing_point
(MP mp
, mp_number
*ret
, mp_number aa
, mp_number bb
, mp_number cc
) {
1225 double d
; /* recursive counter
*/
1226 mpfr_t x
, xx
, x0
, x1
, x2
; /* temporary registers for bisection
*/
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
))
1234 if
(!mpfr_negative_p
(c
)) {
1235 if
(!mpfr_negative_p
(b
)) {
1236 if
(mpfr_positive_p
(c
)) {
1238 } else if
(mpfr_zero_p
(a
) && mpfr_zero_p(b)) {
1246 } else if
(mpfr_zero_p
(a
)) {
1247 if
(!mpfr_positive_p
(b
))
1251 /* Use bisection to find the crossing point...
*/
1253 mpfr_set
(x0
, a
, ROUNDING
);
1254 mpfr_sub
(x1
,a
, b
, ROUNDING
);
1255 mpfr_sub
(x2
,b
, c
, ROUNDING
);
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
);
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
);
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
))
1279 mpfr_set
(x1
,x
, ROUNDING
);
1280 d
= d
+ d
+ epsilonf
;
1283 } while
(d
< fraction_one
);
1284 mpfr_set_d
(scratch
, d
, ROUNDING
);
1285 mpfr_sub
(ret-
>data.num
,scratch
, fraction_one_mpfr_t
, ROUNDING
);
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
));
1291 mpfr_clears
(a
,b
,c
, x
,xx
,x0
,x1
,x2
, scratch
, (mpfr_ptr
)0);
1292 mp_check_mpfr_t
(mp
, ret-
>data.num
);
1297 @ We conclude this set of elementary routines with some simple rounding
1298 and truncation operations.
1301 @ |round_unscaled| rounds a |scaled| and converts it to |int|
1303 int mp_round_unscaled
(mp_number x_orig
) {
1304 double xx
= mp_number_to_double
(x_orig
);
1305 int x
= (int
)ROUND(xx
);
1309 @ |number_floor| floors a number
1312 void mp_number_floor
(mp_number
*i
) {
1313 mpfr_rint_floor
(i-
>data.num
, i-
>data.num
, MPFR_RNDD
);
1316 @ |fraction_to_scaled| rounds a |fraction| and converts it to |scaled|
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
);
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.
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@
>;
1336 mpfr_sqrt
(ret-
>data.num
, x_orig.data.num
, ROUNDING
);
1338 mp_check_mpfr_t
(mp
, ret-
>data.num
);
1342 @ @
<Handle square root of zero...@
>=
1344 if
(mpfr_negative_p
((mpfr_ptr
)x_orig.data.num
)) {
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.",
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
);
1353 @.Square root...replaced by
0@
>;
1354 mp_error
(mp
, msg
, hlp
, true
);
1356 mpfr_set_zero
(ret-
>data.num
,1); /* 1 == positive
*/
1361 @ Pythagorean addition $\psqrt
{a^
2+b^
2}$ is implemented by a quick hack
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);
1377 @ Here is a similar algorithm for $\psqrt
{a^
2-b^
2}$. Same quick hack
, also.
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|@
>;
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
);
1393 mpfr_set
(ret-
>data.num
, a
, ROUNDING
);
1394 mp_check_mpfr_t
(mp
, ret-
>data.num
);
1398 @ @
<Handle erroneous |pyth_sub| and set |a
:=0|@
>=
1400 if
(mpfr_less_p
(a
, b
)) {
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.",
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
);
1412 mp_error
(mp
, msg
, hlp
, true
);
1414 mpfr_set_zero
(a
,1); /* 1 == positive
*/
1418 @ Here is the routine that calculates $
2^
8$ times the natural logarithm
1419 of a |scaled| quantity
;
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@
>;
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
);
1430 mp_check_mpfr_t
(mp
, ret-
>data.num
);
1433 @ @
<Handle non-positive logarithm@
>=
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.",
1440 char
*xstr
= mp_binary_number_tostring
(mp
, x_orig
);
1441 mp_snprintf
(msg
, 256, "Logarithm of %s has been replaced by 0", xstr
);
1443 @.Logarithm...replaced by
0@
>;
1444 mp_error
(mp
, msg
, hlp
, true
);
1445 mpfr_set_zero
(ret-
>data.num
,1); /* 1 == positive
*/
1449 @ Conversely
, the exponential routine calculates $\exp
(x
/2^
8)$
,
1450 when |x| is |scaled|.
1453 void mp_binary_m_exp
(MP mp
, mp_number
*ret
, mp_number x_orig
) {
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
);
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
)$.
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@
>;
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
);
1485 mp_check_mpfr_t
(mp
, ret-
>data.num
);
1489 @ @
<Handle undefined arg@
>=
1491 const char
*hlp
[] = {
1492 "The `angle' between two identical points is undefined.",
1493 "I'm zeroing this one. Proceed, with fingers crossed.",
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
*/
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|.
1505 @ Calculate sines and cosines.
1508 void mp_binary_sin_cos
(MP mp
, mp_number z_orig
, mp_number
*n_cos
, mp_number
*n_sin
) {
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
);
1517 mpfr_sin
(n_sin-
>data.num
, rad
, ROUNDING
);
1518 mpfr_cos
(n_cos-
>data.num
, rad
, ROUNDING
);
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
);
1525 mpfr_clear
(one_eighty
);
1528 @ This is the http
://www-cs-faculty.stanford.edu
/~uno
/programs
/rng.c
1529 with small cosmetic modifications.
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 */
1537 static long ran_x
[KK
]; /* the generator state
*/
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
) */
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
]);
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()" */
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 */
1558 #define TT
70 /* guaranteed separation between streams
*/
1559 #define is_odd
(x
) ((x
)&1) /* units bit of x */
1561 static void ran_start
(long seed
) /* do this before using ran_array
*/
1562 /* long seed selector for different streams
*/
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
*/
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
]);
1582 if
(ss
) ss
>>=1; else t--
;
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;
1590 #define ran_arr_next
() (*ran_arr_ptr
>=0?
*ran_arr_ptr
++: ran_arr_cycle
())
1591 static long ran_arr_cycle
(void
)
1593 if
(ran_arr_ptr
==&ran_arr_dummy)
1594 ran_start
(314159L); /* the user forgot to initialize
*/
1595 ran_array
(ran_arr_buf
,QUALITY
);
1597 ran_arr_ptr
=ran_arr_buf
+1;
1598 return ran_arr_buf
[0];
1604 @ To initialize the |randoms| table
, we call the following routine.
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|
*/
1611 while
(j
>= fraction_one
) {
1615 for
(i
= 0; i
<= 54; i
++) {
1621 mpfr_set_si
(mp-
>randoms
[(i
* 21) % 55].data.num
, j
, ROUNDING
);
1623 mp_new_randoms
(mp
);
1624 mp_new_randoms
(mp
);
1625 mp_new_randoms
(mp
); /* ``warm up'' the array
*/
1627 ran_start
((unsigned long
)seed
);
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
);
1636 @ To consume a random integer for the uniform generator
, the program below will say `|next_unif_random|'.
1639 static void mp_next_unif_random
(MP mp
, mp_number
*ret
) {
1641 unsigned long int op
;
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
);
1654 @ To consume a random fraction
, the program below will say `|next_random|'.
1657 static void mp_next_random
(MP mp
, mp_number
*ret
) {
1658 if
( mp-
>j_random
==0 )
1661 mp-
>j_random
= mp-
>j_random-1
;
1662 mp_number_clone
(ret
, mp-
>randoms
[mp-
>j_random
]);
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.
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.
1673 static void mp_binary_m_unif_rand
(MP mp
, mp_number
*ret
, mp_number x_orig
) {
1674 mp_number y
; /* trial value
*/
1677 char
*r
;mpfr_exp_t e
;
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
);
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
);
1693 mp_number_clone
(ret
, y
);
1694 mp_number_negate
(ret
);
1696 r
= mpfr_get_str
(NULL, // char
*str
,
1697 &e, // mpfr_exp_t *expptr,
1700 ret-
>data.num
, // mpfr_t op
,
1701 ROUNDING
// mpfr_rnd_t rnd
1703 printf
("\nret=%s e=%ld\n",r
,e
);
1705 free_number
(abs_x
);
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\
/}).
1717 static void mp_binary_m_norm_rand
(MP mp
, mp_number
*ret
) {
1723 new_number
(ab_vs_cd
);
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);
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
);
1752 free_number
(abs_x
);
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.
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
) {
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
);
1777 mpfr_mul
(ab
,a
,b
, ROUNDING
);
1778 mpfr_mul
(cd
,c
,d
, ROUNDING
);
1780 mpfr_set
(ret-
>data.num
, zero
, ROUNDING
);
1781 cmp
= mpfr_cmp
(ab
,cd
);
1784 mpfr_set
(ret-
>data.num
, one
, ROUNDING
);
1786 mpfr_set
(ret-
>data.num
, minusone
, ROUNDING
);
1788 mp_check_mpfr_t
(mp
, ret-
>data.num
);
1789 mpfr_clears
(a
,b
,c
,d
,ab
,cd
,(mpfr_ptr
)0);