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
*/
31 @
<Internal library declarations@
>;
34 @
* Math initialization.
36 First
, here are some very important constants.
39 @d E_STRING
"2.7182818284590452353602874713526624977572470936999595749669676277240766303535"
40 @d PI_STRING
"3.1415926535897932384626433832795028841971693993751058209749445923078164062862"
41 @d fraction_multiplier
4096
42 @d angle_multiplier
16
44 @ Here are the functions that are static as they are not used elsewhere
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
);
122 @ We do not want special numbers as return values for functions
, so
:
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);
131 int binary_number_check
(mpfr_t dec
)
134 if
(!mpfr_number_p
(dec
)) {
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
);
142 mpfr_set_zero
(dec
,1); /* 1 == positive
*/
148 void mp_check_mpfr_t
(MP mp
, mpfr_t dec
)
150 mp-
>arith_error
= binary_number_check
(dec
);
156 @ Precision IO uses |double| because |MPFR_PREC_MAX| overflows int.
159 static double precision_bits
;
160 mpfr_prec_t precision_digits_to_bits
(double i
)
164 double precision_bits_to_digits
(mpfr_prec_t d
)
170 @ And these are the ones that
{\it are
} used elsewhere
172 @
<Internal library declarations@
>=
173 void
* mp_initialize_binary_math
(MP mp
);
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
*/
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)
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
;
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
);
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);
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.
243 @d MAX_PRECISION
1000.0
244 @d DEF_PRECISION
34.0
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
();
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);
283 mp_new_number
(mp
, &math->arc_tol_k, mp_fraction_type);
285 mpfr_div_si
(math-
>arc_tol_k.data.num
, one
, 4096, ROUNDING
);
286 /* quit when change in arc length estimate reaches this
*/
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
);
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);
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$
*/
309 mp_new_number
(mp
, &math->twelve_ln_2_k, mp_fraction_type);
311 mpfr_set_d
(math-
>twelve_ln_2_k.data.num
, 139548959.6165 / 65536.0, ROUNDING
);
312 /* $
2^
{24}\cdot12\ln2\approx139548959.6165$
*/
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);
320 mpfr_set_d
(math-
>twelvebits_3.data.num
, 1365 / 65536.0, ROUNDING
);
321 /* $
1365\approx
2^
{12}/3$
*/
323 mp_new_number
(mp
, &math->twentysixbits_sqrt2_t, mp_fraction_type);
325 mpfr_set_d
(math-
>twentysixbits_sqrt2_t.data.num
, 94906265.62 / 65536.0, ROUNDING
);
326 /* $
2^
{26}\sqrt2\approx94906265.62$
*/
328 mp_new_number
(mp
, &math->twentyeightbits_d_t, mp_fraction_type);
330 mpfr_set_d
(math-
>twentyeightbits_d_t.data.num
, 35596754.69 / 65536.0, ROUNDING
);
331 /* $
2^
{28}d\approx35596754.69$
*/
333 mp_new_number
(mp
, &math->twentysevenbits_sqrt2_d_t, mp_fraction_type);
335 mpfr_set_d
(math-
>twentysevenbits_sqrt2_d_t.data.num
, 25170706.63 / 65536.0, ROUNDING
);
336 /* $
2^
{27}\sqrt2\
,d\approx25170706.63$
*/
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
);
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
;
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
);
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
();
457 @ Creating an destroying |mp_number| objects
460 void mp_new_number
(MP mp
, mp_number
*n
, mp_number_type t
) {
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
*/
471 void mp_free_number
(MP mp
, mp_number
*n
) {
474 mpfr_clear
(n-
>data.num
);
477 n-
>type
= mp_nan_type
;
480 @ Here are the low-level functions on |mp_number| items
, setters first.
483 void mp_set_binary_from_int
(mp_number
*A
, int B
) {
484 mpfr_set_si
(A-
>data.num
,B
, ROUNDING
);
486 void mp_set_binary_from_boolean
(mp_number
*A
, int B
) {
487 mpfr_set_si
(A-
>data.num
,B
, ROUNDING
);
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
);
493 void mp_set_binary_from_double
(mp_number
*A
, double B
) {
494 mpfr_set_d
(A-
>data.num
, B
, ROUNDING
);
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
);
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
);
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
);
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
);
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
);
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
);
514 void mp_set_binary_from_of_the_way
(MP mp
, mp_number
*A
, mp_number t
, mp_number B
, mp_number C
) {
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
);
523 mp_check_mpfr_t
(mp
, A-
>data.num
);
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
);
529 void mp_number_add
(mp_number
*A
, mp_number B
) {
530 mpfr_add
(A-
>data.num
,A-
>data.num
,B.data.num
, ROUNDING
);
532 void mp_number_substract
(mp_number
*A
, mp_number B
) {
533 mpfr_sub
(A-
>data.num
,A-
>data.num
,B.data.num
, ROUNDING
);
535 void mp_number_half
(mp_number
*A
) {
536 mpfr_div_si
(A-
>data.num
, A-
>data.num
, 2, ROUNDING
);
538 void mp_number_halfp
(mp_number
*A
) {
539 mpfr_div_si
(A-
>data.num
,A-
>data.num
, 2, ROUNDING
);
541 void mp_number_double
(mp_number
*A
) {
542 mpfr_mul_si
(A-
>data.num
,A-
>data.num
, 2, ROUNDING
);
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
);
547 void mp_number_multiply_int
(mp_number
*A
, int B
) {
548 mpfr_mul_si
(A-
>data.num
,A-
>data.num
, B
, ROUNDING
);
550 void mp_number_divide_int
(mp_number
*A
, int B
) {
551 mpfr_div_si
(A-
>data.num
,A-
>data.num
, B
, ROUNDING
);
553 void mp_binary_abs
(mp_number
*A
) {
554 mpfr_abs
(A-
>data.num
, A-
>data.num
, ROUNDING
);
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
);
560 void mp_number_swap
(mp_number
*A
, mp_number
*B
) {
561 mpfr_swap
(A-
>data.num
, B-
>data.num
);
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
);
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
);
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
);
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
);
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
!
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);
595 @d odd
(A
) (abs
(A
)%2==1)
598 int mp_number_to_int
(mp_number A
) {
600 if
(mpfr_fits_sint_p
(A.data.num
, ROUNDING
)) {
601 result
= mpfr_get_si
(A.data.num
, ROUNDING
);
605 int mp_number_to_boolean
(mp_number A
) {
607 if
(mpfr_fits_sint_p
(A.data.num
, ROUNDING
)) {
608 result
= mpfr_get_si
(A.data.num
, ROUNDING
);
610 return
(result ?
1 : 0);
612 double mp_number_to_double
(mp_number A
) {
614 if
(mpfr_number_p
(A.data.num
)) {
615 res
= mpfr_get_d
(A.data.num
, ROUNDING
);
619 int mp_number_odd
(mp_number A
) {
620 return odd
(mp_number_to_int
(A
));
622 int mp_number_equal
(mp_number A
, mp_number B
) {
623 return mpfr_equal_p
(A.data.num
,B.data.num
);
625 int mp_number_greater
(mp_number A
, mp_number B
) {
626 return mpfr_greater_p
(A.data.num
,B.data.num
);
628 int mp_number_less
(mp_number A
, mp_number B
) {
629 return mpfr_less_p
(A.data.num
,B.data.num
);
631 int mp_number_nonequalabs
(mp_number A
, mp_number B
) {
632 return
!(mpfr_cmpabs
(A.data.num
, B.data.num
)==0);
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.
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|'.
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
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.
658 char
* mp_binnumber_tostring
(mpfr_t n
) {
659 char
*str
= NULL, *buffer
= NULL;
662 if
((str
= mpfr_get_str
(NULL, &exp, 10, 0, n, ROUNDING))>0) {
663 int numprecdigits
= precision_bits_to_digits
(precision_bits
);
667 while
(strlen
(str
)>0 && *(str+strlen(str)-1) == '0' ) {
668 *(str
+strlen
(str
)-1) = '\
0'
; /* get rid of trailing zeroes
*/
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
682 if
(strlen
(str
+j
) == 0) {
686 if
(exp
<=numprecdigits
&& exp > -6) {
688 buffer
[i
++] = str
[j
++];
690 buffer
[i
++] = (str
[j
] ? str
[j
++] : '
0'
);
695 buffer
[i
++] = str
[j
++];
703 while
(absexp--
> 0) {
707 buffer
[i
++] = str
[j
++];
711 buffer
[i
++] = str
[j
++];
715 buffer
[i
++] = str
[j
++];
721 mp_snprintf
(msg
, 256, "%s%d", (exp
>0?
"+":""), (int
)(exp
>0 ?
(exp-1
) : (exp-1
)));
724 buffer
[i
++] = msg
[k
++];
735 char
* mp_binary_number_tostring
(MP mp
, mp_number n
) {
736 return mp_binnumber_tostring
(n.data.num
);
741 void mp_binary_print_number
(MP mp
, mp_number n
) {
742 char
*str
= mp_binary_number_tostring
(mp
, n
);
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
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
);
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.
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.
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.
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@
>
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.
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
);
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
);
808 void mp_binary_make_fraction
(MP mp
, mpfr_t ret
, mpfr_t p
, mpfr_t q
);
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
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@
>
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
);
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
);
830 void mp_binary_take_fraction
(MP mp
, mpfr_t ret
, mpfr_t p
, mpfr_t q
);
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$.
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.
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
);
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.
)
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
);
861 @d halfp
(A
) (integer
)((unsigned
)(A
) >> 1)
863 @
* Scanning numbers in the input
865 The definitions below are temporarily here
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
)
871 static void mp_wrapup_numeric_token
(MP mp
, unsigned char
*start
, unsigned char
*stop
);
873 @ Precision check is TODO
876 void mp_wrapup_numeric_token
(MP mp
, unsigned char
*start
, unsigned char
*stop
) {
879 size_t l
= stop-start
+1;
880 char
*buf
= mp_xmalloc
(mp
, l
+1, 1);
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
));
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
)) {
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.)",
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
);
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:",
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
));
913 set_cur_cmd
((mp_variable_type
)mp_numeric_token
);
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--
;
928 if
(mp-
>buffer
[mp-
>cur_input.loc_field
] == '
+' ||
929 mp-
>buffer
[mp-
>cur_input.loc_field
] == '
-'
) {
930 mp-
>cur_input.loc_field
++;
932 while
(mp-
>char_class
[mp-
>buffer
[mp-
>cur_input.loc_field
]] == digit_class
) {
933 mp-
>cur_input.loc_field
++;
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];
940 while
(mp-
>char_class
[mp-
>buffer
[mp-
>cur_input.loc_field
]] == digit_class
) {
941 mp-
>cur_input.loc_field
++;
944 stop
= &mp->buffer[mp->cur_input.loc_field-1];
945 mp_wrapup_numeric_token
(mp
, start
, stop
);
949 @ We just have to collect bytes.
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];
955 while
(mp-
>char_class
[mp-
>buffer
[mp-
>cur_input.loc_field
]] == digit_class
) {
956 mp-
>cur_input.loc_field
++;
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
++;
966 stop
= &mp->buffer[mp->cur_input.loc_field-1];
967 mp_wrapup_numeric_token
(mp
, start
, stop
);
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
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
)
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.
)
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.
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.
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
*/
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
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
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
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
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
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
1047 if
(!mpfr_equal_p
(t.data.num
, one
)) { // t
!= 1
1048 mpfr_div
(num
, num
, t.data.num
, ROUNDING
); // num
= num
/ t
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
);
1055 mp_binary_make_fraction
(mp
, ret-
>data.num
, num
, denom
);
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
);
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.
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
*/
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|@
>;
1080 mpfr_div
(q
,a
,d
, ROUNDING
);
1081 mpfr_div
(r
,c
,b
, ROUNDING
);
1082 cmp
= mpfr_cmp
(q
,r
);
1085 mpfr_set
(ret-
>data.num
, one
, ROUNDING
);
1087 mpfr_set
(ret-
>data.num
, minusone
, ROUNDING
);
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
);
1097 mpfr_set
(ret-
>data.num
, one
, ROUNDING
);
1101 if
(mpfr_zero_p
(q
)) {
1102 mpfr_set
(ret-
>data.num
, minusone
, ROUNDING
);
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|
*/
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
));
1116 mp_check_mpfr_t
(mp
, ret-
>data.num
);
1117 mpfr_clears
(q
,r
,test
,a
,b
,c
,d
,(mpfr_ptr
)0);
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
);
1127 if
(mpfr_negative_p
(c
)) {
1128 mpfr_neg
(c
, c
, ROUNDING
);
1129 mpfr_neg
(d
, d
, ROUNDING
);
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
);
1136 mpfr_set
(ret-
>data.num
, one
, ROUNDING
);
1139 if
(mpfr_zero_p
(d
)) {
1141 mpfr_set
(ret-
>data.num
, zero
, ROUNDING
);
1143 mpfr_set
(ret-
>data.num
, minusone
, ROUNDING
);
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
);
1158 mpfr_set
(ret-
>data.num
, zero
, ROUNDING
);
1160 mpfr_set
(ret-
>data.num
, minusone
, ROUNDING
);
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.
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$.
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$.
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}$.
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; }
1198 static void mp_binary_crossing_point
(MP mp
, mp_number
*ret
, mp_number aa
, mp_number bb
, mp_number cc
) {
1200 double d
; /* recursive counter
*/
1201 mpfr_t x
, xx
, x0
, x1
, x2
; /* temporary registers for bisection
*/
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
))
1209 if
(!mpfr_negative_p
(c
)) {
1210 if
(!mpfr_negative_p
(b
)) {
1211 if
(mpfr_positive_p
(c
)) {
1213 } else if
(mpfr_zero_p
(a
) && mpfr_zero_p(b)) {
1221 } else if
(mpfr_zero_p
(a
)) {
1222 if
(!mpfr_positive_p
(b
))
1226 /* Use bisection to find the crossing point...
*/
1228 mpfr_set
(x0
, a
, ROUNDING
);
1229 mpfr_sub
(x1
,a
, b
, ROUNDING
);
1230 mpfr_sub
(x2
,b
, c
, ROUNDING
);
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
);
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
);
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
))
1254 mpfr_set
(x1
,x
, ROUNDING
);
1255 d
= d
+ d
+ epsilonf
;
1258 } while
(d
< fraction_one
);
1259 mpfr_set_d
(scratch
, d
, ROUNDING
);
1260 mpfr_sub
(ret-
>data.num
,scratch
, fraction_one_mpfr_t
, ROUNDING
);
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
));
1266 mpfr_clears
(a
,b
,c
, x
,xx
,x0
,x1
,x2
, scratch
, (mpfr_ptr
)0);
1267 mp_check_mpfr_t
(mp
, ret-
>data.num
);
1272 @ We conclude this set of elementary routines with some simple rounding
1273 and truncation operations.
1276 @ |round_unscaled| rounds a |scaled| and converts it to |int|
1278 int mp_round_unscaled
(mp_number x_orig
) {
1279 double xx
= mp_number_to_double
(x_orig
);
1280 int x
= (int
)ROUND(xx
);
1284 @ |number_floor| floors a number
1287 void mp_number_floor
(mp_number
*i
) {
1288 mpfr_rint_floor
(i-
>data.num
, i-
>data.num
, MPFR_RNDD
);
1291 @ |fraction_to_scaled| rounds a |fraction| and converts it to |scaled|
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
);
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.
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@
>;
1311 mpfr_sqrt
(ret-
>data.num
, x_orig.data.num
, ROUNDING
);
1313 mp_check_mpfr_t
(mp
, ret-
>data.num
);
1317 @ @
<Handle square root of zero...@
>=
1319 if
(mpfr_negative_p
((mpfr_ptr
)x_orig.data.num
)) {
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.",
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
);
1328 @.Square root...replaced by
0@
>;
1329 mp_error
(mp
, msg
, hlp
, true
);
1331 mpfr_set_zero
(ret-
>data.num
,1); /* 1 == positive
*/
1336 @ Pythagorean addition $\psqrt
{a^
2+b^
2}$ is implemented by a quick hack
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);
1352 @ Here is a similar algorithm for $\psqrt
{a^
2-b^
2}$. Same quick hack
, also.
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|@
>;
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
);
1368 mpfr_set
(ret-
>data.num
, a
, ROUNDING
);
1369 mp_check_mpfr_t
(mp
, ret-
>data.num
);
1373 @ @
<Handle erroneous |pyth_sub| and set |a
:=0|@
>=
1375 if
(mpfr_less_p
(a
, b
)) {
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.",
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
);
1387 mp_error
(mp
, msg
, hlp
, true
);
1389 mpfr_set_zero
(a
,1); /* 1 == positive
*/
1393 @ Here is the routine that calculates $
2^
8$ times the natural logarithm
1394 of a |scaled| quantity
;
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@
>;
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
);
1405 mp_check_mpfr_t
(mp
, ret-
>data.num
);
1408 @ @
<Handle non-positive logarithm@
>=
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.",
1415 char
*xstr
= mp_binary_number_tostring
(mp
, x_orig
);
1416 mp_snprintf
(msg
, 256, "Logarithm of %s has been replaced by 0", xstr
);
1418 @.Logarithm...replaced by
0@
>;
1419 mp_error
(mp
, msg
, hlp
, true
);
1420 mpfr_set_zero
(ret-
>data.num
,1); /* 1 == positive
*/
1424 @ Conversely
, the exponential routine calculates $\exp
(x
/2^
8)$
,
1425 when |x| is |scaled|.
1428 void mp_binary_m_exp
(MP mp
, mp_number
*ret
, mp_number x_orig
) {
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
);
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
)$.
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@
>;
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
);
1460 mp_check_mpfr_t
(mp
, ret-
>data.num
);
1464 @ @
<Handle undefined arg@
>=
1466 const char
*hlp
[] = {
1467 "The `angle' between two identical points is undefined.",
1468 "I'm zeroing this one. Proceed, with fingers crossed.",
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
*/
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|.
1480 @ Calculate sines and cosines.
1483 void mp_binary_sin_cos
(MP mp
, mp_number z_orig
, mp_number
*n_cos
, mp_number
*n_sin
) {
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
);
1492 mpfr_sin
(n_sin-
>data.num
, rad
, ROUNDING
);
1493 mpfr_cos
(n_cos-
>data.num
, rad
, ROUNDING
);
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
);
1500 mpfr_clear
(one_eighty
);
1503 @ To initialize the |randoms| table
, we call the following routine.
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|
*/
1510 while
(j
>= fraction_one
) {
1514 for
(i
= 0; i
<= 54; i
++) {
1520 mpfr_set_si
(mp-
>randoms
[(i
* 21) % 55].data.num
, j
, ROUNDING
);
1522 mp_new_randoms
(mp
);
1523 mp_new_randoms
(mp
);
1524 mp_new_randoms
(mp
); /* ``warm up'' the array
*/
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
);
1534 @ To consume a random fraction
, the program below will say `|next_random|'.
1537 static void mp_next_random
(MP mp
, mp_number
*ret
) {
1538 if
( mp-
>j_random
==0 )
1541 mp-
>j_random
= mp-
>j_random-1
;
1542 mp_number_clone
(ret
, mp-
>randoms
[mp-
>j_random
]);
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\
/}).
1551 static void mp_binary_m_norm_rand
(MP mp
, mp_number
*ret
) {
1557 new_number
(ab_vs_cd
);
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);
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
);
1586 free_number
(abs_x
);
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.
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
) {
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
);
1611 mpfr_mul
(ab
,a
,b
, ROUNDING
);
1612 mpfr_mul
(cd
,c
,d
, ROUNDING
);
1614 mpfr_set
(ret-
>data.num
, zero
, ROUNDING
);
1615 cmp
= mpfr_cmp
(ab
,cd
);
1618 mpfr_set
(ret-
>data.num
, one
, ROUNDING
);
1620 mpfr_set
(ret-
>data.num
, minusone
, ROUNDING
);
1622 mp_check_mpfr_t
(mp
, ret-
>data.num
);
1623 mpfr_clears
(a
,b
,c
,d
,ab
,cd
,(mpfr_ptr
)0);