1 % $Id
: mpmath.w
2070 2015-10-06 10:35:23Z luigi $
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 % Here is TeX material that gets inserted after \input webmac
9 \font\tenlogo
=logo10
% font used for the METAFONT logo
11 \def\MF
{{\tenlogo META
}\
-{\tenlogo
FONT}}
12 \def\MP
{{\tenlogo META
}\
-{\tenlogo POST
}}
14 \def\title
{Math support functions for
32-bit integer math
}
20 #include
<w2c
/config.h
>
25 #include
"mpmath.h" /* internal header
*/
35 #include
"mpmp.h" /* internal header
*/
36 @
<Internal library declarations@
>;
39 @
* Math initialization.
41 @ Here are the functions that are static as they are not used elsewhere
44 static void mp_scan_fractional_token
(MP mp
, int n
);
45 static void mp_scan_numeric_token
(MP mp
, int n
);
46 static void mp_ab_vs_cd
(MP mp
, mp_number
*ret
, mp_number a
, mp_number b
, mp_number c
, mp_number d
);
47 static void mp_crossing_point
(MP mp
, mp_number
*ret
, mp_number a
, mp_number b
, mp_number c
);
48 static void mp_number_modulo
(mp_number
*a
, mp_number b
);
49 static void mp_print_number
(MP mp
, mp_number n
);
50 static char
* mp_number_tostring
(MP mp
, mp_number n
);
51 static void mp_slow_add
(MP mp
, mp_number
*ret
, mp_number x_orig
, mp_number y_orig
);
52 static void mp_square_rt
(MP mp
, mp_number
*ret
, mp_number x_orig
);
53 static void mp_n_sin_cos
(MP mp
, mp_number z_orig
, mp_number
*n_cos
, mp_number
*n_sin
);
54 static void mp_init_randoms
(MP mp
, int seed
);
55 static void mp_number_angle_to_scaled
(mp_number
*A
);
56 static void mp_number_fraction_to_scaled
(mp_number
*A
);
57 static void mp_number_scaled_to_fraction
(mp_number
*A
);
58 static void mp_number_scaled_to_angle
(mp_number
*A
);
59 static void mp_m_unif_rand
(MP mp
, mp_number
*ret
, mp_number x_orig
);
60 static void mp_m_norm_rand
(MP mp
, mp_number
*ret
);
61 static void mp_m_exp
(MP mp
, mp_number
*ret
, mp_number x_orig
);
62 static void mp_m_log
(MP mp
, mp_number
*ret
, mp_number x_orig
);
63 static void mp_pyth_sub
(MP mp
, mp_number
*r
, mp_number a
, mp_number b
);
64 static void mp_n_arg
(MP mp
, mp_number
*ret
, mp_number x
, mp_number y
);
65 static void mp_velocity
(MP mp
, mp_number
*ret
, mp_number st
, mp_number ct
, mp_number sf
, mp_number cf
, mp_number t
);
66 static void mp_set_number_from_int
(mp_number
*A
, int B
);
67 static void mp_set_number_from_boolean
(mp_number
*A
, int B
);
68 static void mp_set_number_from_scaled
(mp_number
*A
, int B
);
69 static void mp_set_number_from_boolean
(mp_number
*A
, int B
);
70 static void mp_set_number_from_addition
(mp_number
*A
, mp_number B
, mp_number C
);
71 static void mp_set_number_from_substraction
(mp_number
*A
, mp_number B
, mp_number C
);
72 static void mp_set_number_from_div
(mp_number
*A
, mp_number B
, mp_number C
);
73 static void mp_set_number_from_mul
(mp_number
*A
, mp_number B
, mp_number C
);
74 static void mp_set_number_from_int_div
(mp_number
*A
, mp_number B
, int C
);
75 static void mp_set_number_from_int_mul
(mp_number
*A
, mp_number B
, int C
);
76 static void mp_set_number_from_of_the_way
(MP mp
, mp_number
*A
, mp_number t
, mp_number B
, mp_number C
);
77 static void mp_number_negate
(mp_number
*A
);
78 static void mp_number_add
(mp_number
*A
, mp_number B
);
79 static void mp_number_substract
(mp_number
*A
, mp_number B
);
80 static void mp_number_half
(mp_number
*A
);
81 static void mp_number_halfp
(mp_number
*A
);
82 static void mp_number_double
(mp_number
*A
);
83 static void mp_number_add_scaled
(mp_number
*A
, int B
); /* also for negative B
*/
84 static void mp_number_multiply_int
(mp_number
*A
, int B
);
85 static void mp_number_divide_int
(mp_number
*A
, int B
);
86 static void mp_number_abs
(mp_number
*A
);
87 static void mp_number_clone
(mp_number
*A
, mp_number B
);
88 static void mp_number_swap
(mp_number
*A
, mp_number
*B
);
89 static int mp_round_unscaled
(mp_number x_orig
);
90 static int mp_number_to_scaled
(mp_number A
);
91 static int mp_number_to_boolean
(mp_number A
);
92 static int mp_number_to_int
(mp_number A
);
93 static int mp_number_odd
(mp_number A
);
94 static int mp_number_equal
(mp_number A
, mp_number B
);
95 static int mp_number_greater
(mp_number A
, mp_number B
);
96 static int mp_number_less
(mp_number A
, mp_number B
);
97 static int mp_number_nonequalabs
(mp_number A
, mp_number B
);
98 static void mp_number_floor
(mp_number
*i
);
99 static void mp_fraction_to_round_scaled
(mp_number
*x
);
100 static void mp_number_make_scaled
(MP mp
, mp_number
*r
, mp_number p
, mp_number q
);
101 static void mp_number_make_fraction
(MP mp
, mp_number
*r
, mp_number p
, mp_number q
);
102 static void mp_number_take_fraction
(MP mp
, mp_number
*r
, mp_number p
, mp_number q
);
103 static void mp_number_take_scaled
(MP mp
, mp_number
*r
, mp_number p
, mp_number q
);
104 static void mp_new_number
(MP mp
, mp_number
*n
, mp_number_type t
) ;
105 static void mp_free_number
(MP mp
, mp_number
*n
) ;
106 static void mp_free_scaled_math
(MP mp
);
107 static void mp_scaled_set_precision
(MP mp
);
109 @ And these are the ones that
{\it are
} used elsewhere
111 @
<Internal library declarations@
>=
112 void
* mp_initialize_scaled_math
(MP mp
);
113 void mp_set_number_from_double
(mp_number
*A
, double B
);
114 void mp_pyth_add
(MP mp
, mp_number
*r
, mp_number a
, mp_number b
);
115 double mp_number_to_double
(mp_number A
);
119 @d coef_bound
04525252525 /* |fraction| approximation to
7/3 */
120 @d fraction_threshold
2685 /* a |fraction| coefficient less than this is zeroed
*/
121 @d half_fraction_threshold
1342 /* half of |fraction_threshold|
*/
122 @d scaled_threshold
8 /* a |scaled| coefficient less than this is zeroed
*/
123 @d half_scaled_threshold
4 /* half of |scaled_threshold|
*/
124 @d near_zero_angle
26844
125 @d p_over_v_threshold
0x80000
126 @d equation_threshold
64
127 @d tfm_warn_threshold
4096
131 void
* mp_initialize_scaled_math
(MP mp
) {
132 math_data
*math
= (math_data
*)mp_xmalloc
(mp
,1,sizeof
(math_data
));
134 math-
>allocate
= mp_new_number
;
135 math-
>free
= mp_free_number
;
136 mp_new_number
(mp
, &math->precision_default, mp_scaled_type);
137 math-
>precision_default.data.val
= unity
* 10;
138 mp_new_number
(mp
, &math->precision_max, mp_scaled_type);
139 math-
>precision_max.data.val
= unity
* 10;
140 mp_new_number
(mp
, &math->precision_min, mp_scaled_type);
141 math-
>precision_min.data.val
= unity
* 10;
142 /* here are the constants for |scaled| objects
*/
143 mp_new_number
(mp
, &math->epsilon_t, mp_scaled_type);
144 math-
>epsilon_t.data.val
= 1;
145 mp_new_number
(mp
, &math->inf_t, mp_scaled_type);
146 math-
>inf_t.data.val
= EL_GORDO
;
147 mp_new_number
(mp
, &math->warning_limit_t, mp_scaled_type);
148 math-
>warning_limit_t.data.val
= fraction_one
;
149 mp_new_number
(mp
, &math->one_third_inf_t, mp_scaled_type);
150 math-
>one_third_inf_t.data.val
= one_third_EL_GORDO
;
151 mp_new_number
(mp
, &math->unity_t, mp_scaled_type);
152 math-
>unity_t.data.val
= unity
;
153 mp_new_number
(mp
, &math->two_t, mp_scaled_type);
154 math-
>two_t.data.val
= two
;
155 mp_new_number
(mp
, &math->three_t, mp_scaled_type);
156 math-
>three_t.data.val
= three
;
157 mp_new_number
(mp
, &math->half_unit_t, mp_scaled_type);
158 math-
>half_unit_t.data.val
= half_unit
;
159 mp_new_number
(mp
, &math->three_quarter_unit_t, mp_scaled_type);
160 math-
>three_quarter_unit_t.data.val
= three_quarter_unit
;
161 mp_new_number
(mp
, &math->zero_t, mp_scaled_type);
163 mp_new_number
(mp
, &math->arc_tol_k, mp_fraction_type);
164 math-
>arc_tol_k.data.val
= (unity
/4096); /* quit when change in arc length estimate reaches this
*/
165 mp_new_number
(mp
, &math->fraction_one_t, mp_fraction_type);
166 math-
>fraction_one_t.data.val
= fraction_one
;
167 mp_new_number
(mp
, &math->fraction_half_t, mp_fraction_type);
168 math-
>fraction_half_t.data.val
= fraction_half
;
169 mp_new_number
(mp
, &math->fraction_three_t, mp_fraction_type);
170 math-
>fraction_three_t.data.val
= fraction_three
;
171 mp_new_number
(mp
, &math->fraction_four_t, mp_fraction_type);
172 math-
>fraction_four_t.data.val
= fraction_four
;
174 mp_new_number
(mp
, &math->three_sixty_deg_t, mp_angle_type);
175 math-
>three_sixty_deg_t.data.val
= three_sixty_deg
;
176 mp_new_number
(mp
, &math->one_eighty_deg_t, mp_angle_type);
177 math-
>one_eighty_deg_t.data.val
= one_eighty_deg
;
178 /* various approximations
*/
179 mp_new_number
(mp
, &math->one_k, mp_scaled_type);
180 math-
>one_k.data.val
= 1024;
181 mp_new_number
(mp
, &math->sqrt_8_e_k, mp_scaled_type);
182 math-
>sqrt_8_e_k.data.val
= 112429; /* $
2^
{16}\sqrt
{8/e
}\approx
112428.82793$
*/
183 mp_new_number
(mp
, &math->twelve_ln_2_k, mp_fraction_type);
184 math-
>twelve_ln_2_k.data.val
= 139548960; /* $
2^
{24}\cdot12\ln2\approx139548959.6165$
*/
185 mp_new_number
(mp
, &math->coef_bound_k, mp_fraction_type);
186 math-
>coef_bound_k.data.val
= coef_bound
;
187 mp_new_number
(mp
, &math->coef_bound_minus_1, mp_fraction_type);
188 math-
>coef_bound_minus_1.data.val
= coef_bound
- 1;
189 mp_new_number
(mp
, &math->twelvebits_3, mp_scaled_type);
190 math-
>twelvebits_3.data.val
= 1365; /* $
1365\approx
2^
{12}/3$
*/
191 mp_new_number
(mp
, &math->twentysixbits_sqrt2_t, mp_fraction_type);
192 math-
>twentysixbits_sqrt2_t.data.val
= 94906266; /* $
2^
{26}\sqrt2\approx94906265.62$
*/
193 mp_new_number
(mp
, &math->twentyeightbits_d_t, mp_fraction_type);
194 math-
>twentyeightbits_d_t.data.val
= 35596755; /* $
2^
{28}d\approx35596754.69$
*/
195 mp_new_number
(mp
, &math->twentysevenbits_sqrt2_d_t, mp_fraction_type);
196 math-
>twentysevenbits_sqrt2_d_t.data.val
= 25170707; /* $
2^
{27}\sqrt2\
,d\approx25170706.63$
*/
198 mp_new_number
(mp
, &math->fraction_threshold_t, mp_fraction_type);
199 math-
>fraction_threshold_t.data.val
= fraction_threshold
;
200 mp_new_number
(mp
, &math->half_fraction_threshold_t, mp_fraction_type);
201 math-
>half_fraction_threshold_t.data.val
= half_fraction_threshold
;
202 mp_new_number
(mp
, &math->scaled_threshold_t, mp_scaled_type);
203 math-
>scaled_threshold_t.data.val
= scaled_threshold
;
204 mp_new_number
(mp
, &math->half_scaled_threshold_t, mp_scaled_type);
205 math-
>half_scaled_threshold_t.data.val
= half_scaled_threshold
;
206 mp_new_number
(mp
, &math->near_zero_angle_t, mp_angle_type);
207 math-
>near_zero_angle_t.data.val
= near_zero_angle
;
208 mp_new_number
(mp
, &math->p_over_v_threshold_t, mp_fraction_type);
209 math-
>p_over_v_threshold_t.data.val
= p_over_v_threshold
;
210 mp_new_number
(mp
, &math->equation_threshold_t, mp_scaled_type);
211 math-
>equation_threshold_t.data.val
= equation_threshold
;
212 mp_new_number
(mp
, &math->tfm_warn_threshold_t, mp_scaled_type);
213 math-
>tfm_warn_threshold_t.data.val
= tfm_warn_threshold
;
215 math-
>from_int
= mp_set_number_from_int
;
216 math-
>from_boolean
= mp_set_number_from_boolean
;
217 math-
>from_scaled
= mp_set_number_from_scaled
;
218 math-
>from_double
= mp_set_number_from_double
;
219 math-
>from_addition
= mp_set_number_from_addition
;
220 math-
>from_substraction
= mp_set_number_from_substraction
;
221 math-
>from_oftheway
= mp_set_number_from_of_the_way
;
222 math-
>from_div
= mp_set_number_from_div
;
223 math-
>from_mul
= mp_set_number_from_mul
;
224 math-
>from_int_div
= mp_set_number_from_int_div
;
225 math-
>from_int_mul
= mp_set_number_from_int_mul
;
226 math-
>negate
= mp_number_negate
;
227 math-
>add
= mp_number_add
;
228 math-
>substract
= mp_number_substract
;
229 math-
>half
= mp_number_half
;
230 math-
>halfp
= mp_number_halfp
;
231 math-
>do_double
= mp_number_double
;
232 math-
>abs
= mp_number_abs
;
233 math-
>clone
= mp_number_clone
;
234 math-
>swap
= mp_number_swap
;
235 math-
>add_scaled
= mp_number_add_scaled
;
236 math-
>multiply_int
= mp_number_multiply_int
;
237 math-
>divide_int
= mp_number_divide_int
;
238 math-
>to_int
= mp_number_to_int
;
239 math-
>to_boolean
= mp_number_to_boolean
;
240 math-
>to_scaled
= mp_number_to_scaled
;
241 math-
>to_double
= mp_number_to_double
;
242 math-
>odd
= mp_number_odd
;
243 math-
>equal
= mp_number_equal
;
244 math-
>less
= mp_number_less
;
245 math-
>greater
= mp_number_greater
;
246 math-
>nonequalabs
= mp_number_nonequalabs
;
247 math-
>round_unscaled
= mp_round_unscaled
;
248 math-
>floor_scaled
= mp_number_floor
;
249 math-
>fraction_to_round_scaled
= mp_fraction_to_round_scaled
;
250 math-
>make_scaled
= mp_number_make_scaled
;
251 math-
>make_fraction
= mp_number_make_fraction
;
252 math-
>take_fraction
= mp_number_take_fraction
;
253 math-
>take_scaled
= mp_number_take_scaled
;
254 math-
>velocity
= mp_velocity
;
255 math-
>n_arg
= mp_n_arg
;
256 math-
>m_log
= mp_m_log
;
257 math-
>m_exp
= mp_m_exp
;
258 math-
>m_unif_rand
= mp_m_unif_rand
;
259 math-
>m_norm_rand
= mp_m_norm_rand
;
260 math-
>pyth_add
= mp_pyth_add
;
261 math-
>pyth_sub
= mp_pyth_sub
;
262 math-
>fraction_to_scaled
= mp_number_fraction_to_scaled
;
263 math-
>scaled_to_fraction
= mp_number_scaled_to_fraction
;
264 math-
>scaled_to_angle
= mp_number_scaled_to_angle
;
265 math-
>angle_to_scaled
= mp_number_angle_to_scaled
;
266 math-
>init_randoms
= mp_init_randoms
;
267 math-
>sin_cos
= mp_n_sin_cos
;
268 math-
>slow_add
= mp_slow_add
;
269 math-
>sqrt
= mp_square_rt
;
270 math-
>print
= mp_print_number
;
271 math-
>tostring
= mp_number_tostring
;
272 math-
>modulo
= mp_number_modulo
;
273 math-
>ab_vs_cd
= mp_ab_vs_cd
;
274 math-
>crossing_point
= mp_crossing_point
;
275 math-
>scan_numeric
= mp_scan_numeric_token
;
276 math-
>scan_fractional
= mp_scan_fractional_token
;
277 math-
>free_math
= mp_free_scaled_math
;
278 math-
>set_precision
= mp_scaled_set_precision
;
282 void mp_scaled_set_precision
(MP mp
) {
285 void mp_free_scaled_math
(MP mp
) {
286 free_number
(((math_data
*)mp-
>math
)->epsilon_t
);
287 free_number
(((math_data
*)mp-
>math
)->inf_t
);
288 free_number
(((math_data
*)mp-
>math
)->arc_tol_k
);
289 free_number
(((math_data
*)mp-
>math
)->three_sixty_deg_t
);
290 free_number
(((math_data
*)mp-
>math
)->one_eighty_deg_t
);
291 free_number
(((math_data
*)mp-
>math
)->fraction_one_t
);
292 free_number
(((math_data
*)mp-
>math
)->fraction_half_t
);
293 free_number
(((math_data
*)mp-
>math
)->fraction_three_t
);
294 free_number
(((math_data
*)mp-
>math
)->fraction_four_t
);
295 free_number
(((math_data
*)mp-
>math
)->zero_t
);
296 free_number
(((math_data
*)mp-
>math
)->half_unit_t
);
297 free_number
(((math_data
*)mp-
>math
)->three_quarter_unit_t
);
298 free_number
(((math_data
*)mp-
>math
)->unity_t
);
299 free_number
(((math_data
*)mp-
>math
)->two_t
);
300 free_number
(((math_data
*)mp-
>math
)->three_t
);
301 free_number
(((math_data
*)mp-
>math
)->one_third_inf_t
);
302 free_number
(((math_data
*)mp-
>math
)->warning_limit_t
);
303 free_number
(((math_data
*)mp-
>math
)->one_k
);
304 free_number
(((math_data
*)mp-
>math
)->sqrt_8_e_k
);
305 free_number
(((math_data
*)mp-
>math
)->twelve_ln_2_k
);
306 free_number
(((math_data
*)mp-
>math
)->coef_bound_k
);
307 free_number
(((math_data
*)mp-
>math
)->coef_bound_minus_1
);
308 free_number
(((math_data
*)mp-
>math
)->twelvebits_3
);
309 free_number
(((math_data
*)mp-
>math
)->twentysixbits_sqrt2_t
);
310 free_number
(((math_data
*)mp-
>math
)->twentyeightbits_d_t
);
311 free_number
(((math_data
*)mp-
>math
)->twentysevenbits_sqrt2_d_t
);
312 free_number
(((math_data
*)mp-
>math
)->fraction_threshold_t
);
313 free_number
(((math_data
*)mp-
>math
)->half_fraction_threshold_t
);
314 free_number
(((math_data
*)mp-
>math
)->scaled_threshold_t
);
315 free_number
(((math_data
*)mp-
>math
)->half_scaled_threshold_t
);
316 free_number
(((math_data
*)mp-
>math
)->near_zero_angle_t
);
317 free_number
(((math_data
*)mp-
>math
)->p_over_v_threshold_t
);
318 free_number
(((math_data
*)mp-
>math
)->equation_threshold_t
);
319 free_number
(((math_data
*)mp-
>math
)->tfm_warn_threshold_t
);
323 @ Creating an destroying |mp_number| objects
326 void mp_new_number
(MP mp
, mp_number
*n
, mp_number_type t
) {
334 void mp_free_number
(MP mp
, mp_number
*n
) {
336 n-
>type
= mp_nan_type
;
339 @ Here are the low-level functions on |mp_number| items
, setters first.
342 void mp_set_number_from_int
(mp_number
*A
, int B
) {
345 void mp_set_number_from_boolean
(mp_number
*A
, int B
) {
348 void mp_set_number_from_scaled
(mp_number
*A
, int B
) {
351 void mp_set_number_from_double
(mp_number
*A
, double B
) {
352 A-
>data.val
= (int
)(B
*65536.0);
354 void mp_set_number_from_addition
(mp_number
*A
, mp_number B
, mp_number C
) {
355 A-
>data.val
= B.data.val
+C.data.val
;
357 void mp_set_number_from_substraction
(mp_number
*A
, mp_number B
, mp_number C
) {
358 A-
>data.val
= B.data.val-C.data.val
;
360 void mp_set_number_from_div
(mp_number
*A
, mp_number B
, mp_number C
) {
361 A-
>data.val
= B.data.val
/ C.data.val
;
363 void mp_set_number_from_mul
(mp_number
*A
, mp_number B
, mp_number C
) {
364 A-
>data.val
= B.data.val
* C.data.val
;
366 void mp_set_number_from_int_div
(mp_number
*A
, mp_number B
, int C
) {
367 A-
>data.val
= B.data.val
/ C
;
369 void mp_set_number_from_int_mul
(mp_number
*A
, mp_number B
, int C
) {
370 A-
>data.val
= B.data.val
* C
;
372 void mp_set_number_from_of_the_way
(MP mp
, mp_number
*A
, mp_number t
, mp_number B
, mp_number C
) {
373 A-
>data.val
= B.data.val
- mp_take_fraction
(mp
, (B.data.val
- C.data.val
), t.data.val
);
375 void mp_number_negate
(mp_number
*A
) {
376 A-
>data.val
= -A-
>data.val
;
378 void mp_number_add
(mp_number
*A
, mp_number B
) {
379 A-
>data.val
= A-
>data.val
+ B.data.val
;
381 void mp_number_substract
(mp_number
*A
, mp_number B
) {
382 A-
>data.val
= A-
>data.val
- B.data.val
;
384 void mp_number_half
(mp_number
*A
) {
385 A-
>data.val
= A-
>data.val
/2;
387 void mp_number_halfp
(mp_number
*A
) {
388 A-
>data.val
= (A-
>data.val
>>1);
390 void mp_number_double
(mp_number
*A
) {
391 A-
>data.val
= A-
>data.val
+ A-
>data.val
;
393 void mp_number_add_scaled
(mp_number
*A
, int B
) { /* also for negative B
*/
394 A-
>data.val
= A-
>data.val
+ B
;
396 void mp_number_multiply_int
(mp_number
*A
, int B
) {
397 A-
>data.val
= B
* A-
>data.val
;
399 void mp_number_divide_int
(mp_number
*A
, int B
) {
400 A-
>data.val
= A-
>data.val
/ B
;
402 void mp_number_abs
(mp_number
*A
) {
403 A-
>data.val
= abs
(A-
>data.val
);
405 void mp_number_clone
(mp_number
*A
, mp_number B
) {
406 A-
>data.val
= B.data.val
;
408 void mp_number_swap
(mp_number
*A
, mp_number
*B
) {
409 int swap_tmp
= A-
>data.val
;
410 A-
>data.val
= B-
>data.val
;
411 B-
>data.val
= swap_tmp
;
413 void mp_number_fraction_to_scaled
(mp_number
*A
) {
414 A-
>type
= mp_scaled_type
;
415 A-
>data.val
= A-
>data.val
/ 4096;
417 void mp_number_angle_to_scaled
(mp_number
*A
) {
418 A-
>type
= mp_scaled_type
;
419 if
(A-
>data.val
>= 0) {
420 A-
>data.val
= (A-
>data.val
+ 8) / 16;
422 A-
>data.val
= -((-A-
>data.val
+ 8) / 16);
425 void mp_number_scaled_to_fraction
(mp_number
*A
) {
426 A-
>type
= mp_fraction_type
;
427 A-
>data.val
= A-
>data.val
* 4096;
429 void mp_number_scaled_to_angle
(mp_number
*A
) {
430 A-
>type
= mp_angle_type
;
431 A-
>data.val
= A-
>data.val
* 16;
438 int mp_number_to_int
(mp_number A
) {
441 int mp_number_to_scaled
(mp_number A
) {
444 int mp_number_to_boolean
(mp_number A
) {
447 double mp_number_to_double
(mp_number A
) {
448 return
(A.data.val
/65536.0);
450 int mp_number_odd
(mp_number A
) {
451 return odd
(A.data.val
);
453 int mp_number_equal
(mp_number A
, mp_number B
) {
454 return
(A.data.val
==B.data.val
);
456 int mp_number_greater
(mp_number A
, mp_number B
) {
457 return
(A.data.val
>B.data.val
);
459 int mp_number_less
(mp_number A
, mp_number B
) {
460 return
(A.data.val
<B.data.val
);
462 int mp_number_nonequalabs
(mp_number A
, mp_number B
) {
463 return
(!(abs
(A.data.val
)==abs
(B.data.val
)));
466 @ Fixed-point arithmetic is done on
{\sl scaled integers\
/} that are multiples
467 of $
2^
{-16}$. In other words
, a binary point is assumed to be sixteen bit
468 positions from the right end of a binary computer word.
470 @d unity
0x10000 /* $
2^
{16}$
, represents
1.00000 */
471 @d two
(2*unity
) /* $
2^
{17}$
, represents
2.00000 */
472 @d three
(3*unity
) /* $
2^
{17}+2^
{16}$
, represents
3.00000 */
473 @d half_unit
(unity
/2) /* $
2^
{15}$
, represents
0.50000 */
474 @d three_quarter_unit
(3*(unity
/4)) /* $
3\cdot2^
{14}$
, represents
0.75000 */
476 @d EL_GORDO
0x7fffffff /* $
2^
{31}-1$
, the largest value that \MP\ likes
*/
477 @d one_third_EL_GORDO
05252525252
479 @ One of \MP's most common operations is the calculation of
480 $\lfloor
{a
+b\over2
}\rfloor$
,
481 the midpoint of two given integers |a| and~|b|. The most decent way to do
482 this is to write `|
(a
+b
)/2|'
; but on many machines it is more efficient
483 to calculate `|
(a
+b
)>>1|'.
485 Therefore the midpoint operation will always be denoted by `|half
(a
+b
)|'
486 in this program. If \MP\ is being implemented with languages that permit
487 binary shifting
, the |half| macro should be changed to make this operation
488 as efficient as possible. Since some systems have shift operators that can
489 only be trusted to work on positive numbers
, there is also a macro |halfp|
490 that is used only when the quantity being halved is known to be positive
493 @d halfp
(A
) (integer
)((unsigned
)(A
) >> 1)
495 @ Here is a procedure analogous to |print_int|. If the output
496 of this procedure is subsequently read by \MP\ and converted by the
497 |round_decimals| routine above
, it turns out that the original value will
498 be reproduced exactly. A decimal point is printed only if the value is
499 not an integer. If there is more than one way to print the result with
500 the optimum number of digits following the decimal point
, the closest
501 possible value is given.
503 The invariant relation in the \
&{repeat} loop is that a sequence of
504 decimal digits yet to be printed will yield the original number if and only if
505 they form a fraction~$f$ in the range $s-\delta\L10\cdot2^
{16}f
<s$.
506 We can stop if and only if $f
=0$ satisfies this condition
; the loop will
507 terminate before $s$ can possibly become zero.
510 static void mp_print_scaled
(MP mp
, int s
); /* scaled
*/
511 static char
*mp_string_scaled
(MP mp
, int s
);
514 static void mp_print_scaled
(MP mp
, int s
) { /* s
=scaled prints scaled real
, rounded to five digits
*/
515 int delta
; /* amount of allowable inaccuracy
, scaled
*/
517 mp_print_char
(mp
, xord
('
-'
));
518 s
= -s
; /* print the sign
, if negative
*/
520 mp_print_int
(mp
, s
/ unity
); /* print the integer part
*/
521 s
= 10 * (s
% unity
) + 5;
524 mp_print_char
(mp
, xord
('.'
));
527 s
= s
+ 0100000 - (delta
/ 2); /* round the final digit
*/
528 mp_print_char
(mp
, xord
('
0'
+ (s
/ unity
)));
529 s
= 10 * (s
% unity
);
535 static char
*mp_string_scaled
(MP mp
, int s
) { /* s
=scaled prints scaled real
, rounded to five digits
*/
536 static char scaled_string
[32];
537 int delta
; /* amount of allowable inaccuracy
, scaled
*/
540 scaled_string
[i
++] = xord
('
-'
);
541 s
= -s
; /* print the sign
, if negative
*/
543 /* print the integer part
*/
544 mp_snprintf
((scaled_string
+i
), 12, "%d", (int
) (s
/ unity
));
545 while
(*(scaled_string
+i
)) i
++;
547 s
= 10 * (s
% unity
) + 5;
550 scaled_string
[i
++] = xord
('.'
);
553 s
= s
+ 0100000 - (delta
/ 2); /* round the final digit
*/
554 scaled_string
[i
++] = xord
('
0'
+ (s
/ unity
));
555 s
= 10 * (s
% unity
);
559 scaled_string
[i
] = '\
0'
;
560 return scaled_string
;
563 @ Addition is not always checked to make sure that it doesn't overflow
,
564 but in places where overflow isn't too unlikely the |slow_add| routine
568 void mp_slow_add
(MP mp
, mp_number
*ret
, mp_number x_orig
, mp_number y_orig
) {
573 if
(y
<= EL_GORDO
- x
) {
574 ret-
>data.val
= x
+ y
;
576 mp-
>arith_error
= true
;
577 ret-
>data.val
= EL_GORDO
;
579 } else if
(-y
<= EL_GORDO
+ x
) {
580 ret-
>data.val
= x
+ y
;
582 mp-
>arith_error
= true
;
583 ret-
>data.val
= -EL_GORDO
;
587 @ The |make_fraction| routine produces the |fraction| equivalent of
588 |p
/q|
, given integers |p| and~|q|
; it computes the integer
589 $f
=\lfloor2^
{28}p
/q
+{1\over2
}\rfloor$
, when $p$ and $q$ are
590 positive. If |p| and |q| are both of the same scaled type |t|
,
591 the ``type relation'' |make_fraction
(t
,t
)=fraction| is valid
;
592 and it's also possible to use the subroutine ``backwards
,'' using
593 the relation |make_fraction
(t
,fraction
)=t| between scaled types.
595 If the result would have magnitude $
2^
{31}$ or more
, |make_fraction|
596 sets |arith_error
:=true|. Most of \MP's internal computations have
597 been designed to avoid this sort of error.
599 If this subroutine were programmed in assembly language on a typical
600 machine
, we could simply compute |
(@t$
2^
{28}$@
>*p
)div q|
, since a
601 double-precision product can often be input to a fixed-point division
602 instruction. But when we are restricted to int-eger arithmetic it
603 is necessary either to resort to multiple-precision maneuvering
604 or to use a simple but slow iteration. The multiple-precision technique
605 would be about three times faster than the code adopted here
, but it
606 would be comparatively long and tricky
, involving about sixteen
607 additional multiplications and divisions.
609 This operation is part of \MP's ``inner loop''
; indeed
, it will
610 consume nearly
10\pct
! of the running time
(exclusive of input and output
)
611 if the code below is left unchanged. A machine-dependent recoding
612 will therefore make \MP\ run faster. The present implementation
613 is highly portable
, but slow
; it avoids multiplication and division
614 except in the initial stage. System wizards should be careful to
615 replace it with a routine that is guaranteed to produce identical
616 results in all cases.
617 @^system dependencies@
>
619 As noted below
, a few more routines should also be replaced by machine-dependent
620 code
, for efficiency. But when a procedure is not part of the ``inner loop
,''
621 such changes aren't advisable
; simplicity and robustness are
622 preferable to trickery
, unless the cost is too high.
625 @ We need these preprocessor values
627 @d TWEXP31
2147483648.0
628 @d TWEXP28
268435456.0
630 @d TWEXP_16
(1.0/65536.0)
631 @d TWEXP_28
(1.0/268435456.0)
635 static integer mp_make_fraction
(MP mp
, integer p
, integer q
) {
638 mp_confusion
(mp
, "/");
639 @
:this can't happen
/}{\quad \.
/@
>
642 d
= TWEXP28
* (double
) p
/ (double
) q
;
646 mp-
>arith_error
= true
;
651 if
(d
== (double
) i
&& (((q > 0 ? -q : q) & 077777)
652 * (((i
& 037777) << 1) - 1) & 04000) != 0)
657 mp-
>arith_error
= true
;
662 if
(d
== (double
) i
&& (((q > 0 ? q : -q) & 077777)
663 * (((i
& 037777) << 1) + 1) & 04000) != 0)
670 void mp_number_make_fraction
(MP mp
, mp_number
*ret
, mp_number p
, mp_number q
) {
671 ret-
>data.val
= mp_make_fraction
(mp
, p.data.val
, q.data.val
);
675 @ The dual of |make_fraction| is |take_fraction|
, which multiplies a
676 given integer~|q| by a fraction~|f|. When the operands are positive
, it
677 computes $p
=\lfloor qf
/2^
{28}+{1\over2
}\rfloor$
, a symmetric function
680 This routine is even more ``inner loopy'' than |make_fraction|
;
681 the present implementation consumes almost
20\pct
! of \MP's computation
682 time during typical jobs
, so a machine-language substitute is advisable.
683 @^inner loop@
> @^system dependencies@
>
685 @
<Internal library declarations@
>=
686 /* still in use by tfmin.w
*/
687 integer mp_take_fraction
(MP mp
, integer q
, int f
);
690 integer mp_take_fraction
(MP mp
, integer p
, int q
) { /* q
= fraction
*/
693 d
= (double
) p
*(double
) q
*TWEXP_28
;
697 if
(d
!= TWEXP31 ||
(((p
& 077777) * (q & 077777)) & 040000) == 0)
698 mp-
>arith_error
= true
;
702 if
(d
== (double
) i
&& (((p & 077777) * (q & 077777)) & 040000) != 0)
707 if
(d
!= -TWEXP31 ||
((-(p
& 077777) * (q & 077777)) & 040000) == 0)
708 mp-
>arith_error
= true
;
712 if
(d
== (double
) i
&& ((-(p & 077777) * (q & 077777)) & 040000) != 0)
717 void mp_number_take_fraction
(MP mp
, mp_number
*ret
, mp_number p_orig
, mp_number q_orig
) {
718 ret-
>data.val
= mp_take_fraction
(mp
, p_orig.data.val
, q_orig.data.val
);
722 @ When we want to multiply something by a |scaled| quantity
, we use a scheme
723 analogous to |take_fraction| but with a different scaling.
724 Given positive operands
, |take_scaled|
725 computes the quantity $p
=\lfloor qf
/2^
{16}+{1\over2
}\rfloor$.
727 Once again it is a good idea to use a machine-language replacement if
728 possible
; otherwise |take_scaled| will use more than
2\pct
! of the running time
729 when the Computer Modern fonts are being generated.
733 static integer mp_take_scaled
(MP mp
, integer q
, int f
);
736 static integer mp_take_scaled
(MP mp
, integer p
, int q
) { /* q
= scaled
*/
739 d
= (double
) p
*(double
) q
*TWEXP_16
;
743 if
(d
!= TWEXP31 ||
(((p
& 077777) * (q & 077777)) & 040000) == 0)
744 mp-
>arith_error
= true
;
748 if
(d
== (double
) i
&& (((p & 077777) * (q & 077777)) & 040000) != 0)
753 if
(d
!= -TWEXP31 ||
((-(p
& 077777) * (q & 077777)) & 040000) == 0)
754 mp-
>arith_error
= true
;
758 if
(d
== (double
) i
&& ((-(p & 077777) * (q & 077777)) & 040000) != 0)
763 void mp_number_take_scaled
(MP mp
, mp_number
*ret
, mp_number p_orig
, mp_number q_orig
) {
764 ret-
>data.val
= mp_take_scaled
(mp
, p_orig.data.val
, q_orig.data.val
);
768 @ For completeness
, there's also |make_scaled|
, which computes a
769 quotient as a |scaled| number instead of as a |fraction|.
770 In other words
, the result is $\lfloor2^
{16}p
/q
+{1\over2
}\rfloor$
, if the
771 operands are positive. \
(This procedure is not used especially often
,
772 so it is not part of \MP's inner loop.
)
774 @
<Internal library ...@
>=
775 /* still in use by svgout.w
*/
776 int mp_make_scaled
(MP mp
, integer p
, integer q
);
779 int mp_make_scaled
(MP mp
, integer p
, integer q
) { /* return scaled
*/
782 mp_confusion
(mp
, "/");
783 @
:this can't happen
/}{\quad \.
/@
> {
785 d
= TWEXP16
* (double
) p
/ (double
) q
;
789 mp-
>arith_error
= true
;
793 if
(d
== (double
) i
&& (((q > 0 ? -q : q) & 077777)
794 * (((i
& 037777) << 1) - 1) & 04000) != 0)
799 mp-
>arith_error
= true
;
803 if
(d
== (double
) i
&& (((q > 0 ? q : -q) & 077777)
804 * (((i
& 037777) << 1) + 1) & 04000) != 0)
810 void mp_number_make_scaled
(MP mp
, mp_number
*ret
, mp_number p_orig
, mp_number q_orig
) {
811 ret-
>data.val
= mp_make_scaled
(mp
, p_orig.data.val
, q_orig.data.val
);
814 @ The following function is used to create a scaled integer from a given decimal
815 fraction $
(.d_0d_1\ldots d_
{k-1
})$
, where |
0<=k
<=17|.
818 static int mp_round_decimals
(MP mp
, unsigned char
*b
, quarterword k
);
821 static int mp_round_decimals
(MP mp
, unsigned char
*b
, quarterword k
) { /* return
: scaled
*/
822 /* converts a decimal fraction
*/
823 unsigned a
= 0; /* the accumulator
*/
825 (void
)mp
; /* Will be needed later
*/
826 for
( l
= k-1
; l
>= 0; l--
) {
827 if
(l
<16) /* digits for |k
>=17| cannot affect the result
*/
828 a
= (a
+ (unsigned
) (*(b
+l
) - '
0'
) * two
) / 10;
830 return
(int
) halfp
(a
+ 1);
833 @
* Scanning numbers in the input
835 The definitions below are temporarily here
837 @d set_cur_cmd
(A
) mp-
>cur_mod_-
>type
=(A
)
838 @d set_cur_mod
(A
) mp-
>cur_mod_-
>data.n.data.val
=(A
)
841 static void mp_wrapup_numeric_token
(MP mp
, int n
, int f
);
844 static void mp_wrapup_numeric_token
(MP mp
, int n
, int f
) { /* n
,f
: scaled
*/
845 int mod
; /* scaled
*/
847 mod
= (n
* unity
+ f
);
849 if
(mod
>= fraction_one
) {
850 if
(internal_value
(mp_warning_check
).data.val
> 0 &&
851 (mp-
>scanner_status
!= tex_flushing
)) {
853 const char
*hlp
[] = {"It is at least 4096. Continue and I'll try to cope",
854 "with that big value; but it might be dangerous.",
855 "(Set warningcheck:=0 to suppress this message.)",
857 mp_snprintf
(msg
, 256, "Number is too large (%s)", mp_string_scaled
(mp
,mod
));
858 @.Number is too large@
>;
859 mp_error
(mp
, msg
, hlp
, true
);
862 } else if
(mp-
>scanner_status
!= tex_flushing
) {
863 const char
*hlp
[] = {"I can\'t handle numbers bigger than 32767.99998;",
864 "so I've changed your constant to that maximum amount.",
866 mp_error
(mp
, "Enormous number has been reduced", hlp
, false
);
867 @.Enormous number...@
>;
868 set_cur_mod
(EL_GORDO
);
870 set_cur_cmd
((mp_variable_type
)mp_numeric_token
);
874 void mp_scan_fractional_token
(MP mp
, int n
) { /* n
: scaled
*/
879 mp-
>cur_input.loc_field
++;
880 } while
(mp-
>char_class
[mp-
>buffer
[mp-
>cur_input.loc_field
]] == digit_class
);
881 f
= mp_round_decimals
(mp
, (unsigned char
*)(mp-
>buffer
+mp-
>cur_input.loc_field-k
), (quarterword
) k
);
886 mp_wrapup_numeric_token
(mp
, n
, f
);
891 void mp_scan_numeric_token
(MP mp
, int n
) { /* n
: scaled
*/
892 while
(mp-
>char_class
[mp-
>buffer
[mp-
>cur_input.loc_field
]] == digit_class
) {
894 n
= 10 * n
+ mp-
>buffer
[mp-
>cur_input.loc_field
] - '
0'
;
895 mp-
>cur_input.loc_field
++;
897 if
(!(mp-
>buffer
[mp-
>cur_input.loc_field
] == '.'
&&
898 mp-
>char_class
[mp-
>buffer
[mp-
>cur_input.loc_field
+ 1]] == digit_class
)) {
899 mp_wrapup_numeric_token
(mp
, n
, 0);
901 mp-
>cur_input.loc_field
++;
902 mp_scan_fractional_token
(mp
, n
);
906 @ The |scaled| quantities in \MP\ programs are generally supposed to be
907 less than $
2^
{12}$ in absolute value
, so \MP\ does much of its internal
908 arithmetic with
28~significant bits of precision. A |fraction| denotes
909 a scaled integer whose binary point is assumed to be
28 bit positions
912 @d fraction_half
01000000000 /* $
2^
{27}$
, represents
0.50000000 */
913 @d fraction_one
02000000000 /* $
2^
{28}$
, represents
1.00000000 */
914 @d fraction_two
04000000000 /* $
2^
{29}$
, represents
2.00000000 */
915 @d fraction_three
06000000000 /* $
3\cdot2^
{28}$
, represents
3.00000000 */
916 @d fraction_four
010000000000 /* $
2^
{30}$
, represents
4.00000000 */
918 @ Here is a typical example of how the routines above can be used.
919 It computes the function
920 $$
{1\over3\tau
}f
(\theta
,\phi
)=
921 {\tau^
{-1}\bigl
(2+\sqrt2\
,(\sin\theta-
{1\over16
}\sin\phi
)
922 (\sin\phi-
{1\over16
}\sin\theta
)(\cos\theta-\cos\phi
)\bigr
)\over
923 3\
,\bigl
(1+{1\over2
}(\sqrt5-1
)\cos\theta
+{1\over2
}(3-\sqrt5\
,)\cos\phi\bigr
)},$$
924 where $\tau$ is a |scaled| ``tension'' parameter. This is \MP's magic
925 fudge factor for placing the first control point of a curve that starts
926 at an angle $\theta$ and ends at an angle $\phi$ from the straight path.
927 (Actually
, if the stated quantity exceeds
4, \MP\ reduces it to~
4.
)
929 The trigonometric quantity to be multiplied by $\sqrt2$ is less than $\sqrt2$.
930 (It's a sum of eight terms whose absolute values can be bounded using
931 relations such as $\sin\theta\cos\theta\L
{1\over2
}$.
) Thus the numerator
932 is positive
; and since the tension $\tau$ is constrained to be at least
933 $
3\over4$
, the numerator is less than $
16\over3$. The denominator is
934 nonnegative and at most~
6. Hence the fixed-point calculations below
935 are guaranteed to stay within the bounds of a
32-bit computer word.
937 The angles $\theta$ and $\phi$ are given implicitly in terms of |fraction|
938 arguments |st|
, |ct|
, |sf|
, and |cf|
, representing $\sin\theta$
, $\cos\theta$
,
939 $\sin\phi$
, and $\cos\phi$
, respectively.
942 void mp_velocity
(MP mp
, mp_number
*ret
, mp_number st
, mp_number ct
, mp_number sf
,
943 mp_number cf
, mp_number t
) {
944 integer acc
, num
, denom
; /* registers for intermediate calculations
*/
945 acc
= mp_take_fraction
(mp
, st.data.val
- (sf.data.val
/ 16), sf.data.val
- (st.data.val
/ 16));
946 acc
= mp_take_fraction
(mp
, acc
, ct.data.val
- cf.data.val
);
947 num
= fraction_two
+ mp_take_fraction
(mp
, acc
, 379625062);
948 /* $
2^
{28}\sqrt2\approx379625062.497$
*/
950 fraction_three
+ mp_take_fraction
(mp
, ct.data.val
,
951 497706707) + mp_take_fraction
(mp
, cf.data.val
,
953 /* $
3\cdot2^
{27}\cdot
(\sqrt5-1
)\approx497706706.78$ and
954 $
3\cdot2^
{27}\cdot
(3-\sqrt5\
,)\approx307599661.22$
*/
955 if
(t.data.val
!= unity
)
956 num
= mp_make_scaled
(mp
, num
, t.data.val
); /* |make_scaled
(fraction
,scaled
)=fraction|
*/
957 if
(num
/ 4 >= denom
) {
958 ret-
>data.val
= fraction_four
;
960 ret-
>data.val
= mp_make_fraction
(mp
, num
, denom
);
962 /* printf
("num,denom=%f,%f -=> %f\n", num
/65536.0, denom
/65536.0, ret.data.val
/65536.0);*/
966 @ The following somewhat different subroutine tests rigorously if $ab$ is
967 greater than
, equal to
, or less than~$cd$
,
968 given integers $
(a
,b
,c
,d
)$. In most cases a quick decision is reached.
969 The result is $
+1$
, 0, or~$
-1$ in the three respective cases.
972 static 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
) {
973 integer q
, r
; /* temporary registers
*/
980 @
<Reduce to the case that |a
,c
>=0|
, |b
,d
>0|@
>;
985 ret-
>data.val
= (q
> r ?
1 : -1);
991 ret-
>data.val
= (q ?
1 : 0);
1002 } /* now |a
>d
>0| and |c
>b
>0|
*/
1006 @ @
<Reduce to the case that |a...@
>=
1017 if
((a
== 0 || b
== 0) && (c == 0 || d == 0))
1024 ret-
>data.val
= (a
== 0 ?
0 : -1);
1033 } else if
(b
<= 0) {
1034 if
(b
< 0 && a > 0) {
1038 ret-
>data.val
= (c
== 0 ?
0 : -1);
1042 @ Now here's a subroutine that's handy for all sorts of path computations
:
1043 Given a quadratic polynomial $B
(a
,b
,c
;t
)$
, the |crossing_point| function
1044 returns the unique |fraction| value |t| between
0 and~
1 at which
1045 $B
(a
,b
,c
;t
)$ changes from positive to negative
, or returns
1046 |t
=fraction_one
+1| if no such value exists. If |a
<0|
(so that $B
(a
,b
,c
;t
)$
1047 is already negative at |t
=0|
), |crossing_point| returns the value zero.
1049 The general bisection method is quite simple when $n
=2$
, hence
1050 |crossing_point| does not take much time. At each stage in the
1051 recursion we have a subinterval defined by |l| and~|j| such that
1052 $B
(a
,b
,c
;2^
{-l
}(j
+t
))=B
(x_0
,x_1
,x_2
;t
)$
, and we want to ``zero in'' on
1053 the subinterval where $x_0\G0$ and $\min
(x_1
,x_2
)<0$.
1055 It is convenient for purposes of calculation to combine the values
1056 of |l| and~|j| in a single variable $d
=2^l
+j$
, because the operation
1057 of bisection then corresponds simply to doubling $d$ and possibly
1058 adding~
1. Furthermore it proves to be convenient to modify
1059 our previous conventions for bisection slightly
, maintaining the
1060 variables $X_0
=2^lx_0$
, $X_1
=2^l
(x_0-x_1
)$
, and $X_2
=2^l
(x_1-x_2
)$.
1061 With these variables the conditions $x_0\ge0$ and $\min
(x_1
,x_2
)<0$ are
1062 equivalent to $\max
(X_1
,X_1
+X_2
)>X_0\ge0$.
1064 The following code maintains the invariant relations
1065 $
0\L|x0|
<\max
(|x1|
,|x1|
+|x2|
)$
,
1066 $\vert|x1|\vert
<2^
{30}$
, $\vert|x2|\vert
<2^
{30}$
;
1067 it has been constructed in such a way that no arithmetic overflow
1068 will occur if the inputs satisfy
1069 $a
<2^
{30}$
, $\vert a-b\vert
<2^
{30}$
, and $\vert b-c\vert
<2^
{30}$.
1071 @d no_crossing
{ ret-
>data.val
= fraction_one
+ 1; return
; }
1072 @d one_crossing
{ ret-
>data.val
= fraction_one
; return
; }
1073 @d zero_crossing
{ ret-
>data.val
= 0; return
; }
1076 static void mp_crossing_point
(MP mp
, mp_number
*ret
, mp_number aa
, mp_number bb
, mp_number cc
) {
1078 integer d
; /* recursive counter
*/
1079 integer x
, xx
, x0
, x1
, x2
; /* temporary registers for bisection
*/
1089 } else if
((a
== 0) && (b == 0)) {
1097 } else if
(a
== 0) {
1102 /* Use bisection to find the crossing point...
*/
1129 } while
(d
< fraction_one
);
1130 ret-
>data.val
= (d
- fraction_one
);
1134 @ We conclude this set of elementary routines with some simple rounding
1135 and truncation operations.
1138 @ |round_unscaled| rounds a |scaled| and converts it to |int|
1140 int mp_round_unscaled
(mp_number x_orig
) {
1141 int x
= x_orig.data.val
;
1143 return
1+((x-32768
) / 65536);
1144 } else if
( x
>=-32768) {
1147 return
-(1+((-(x
+1)-32768) / 65536));
1151 @ |number_floor| floors a |scaled|
1154 void mp_number_floor
(mp_number
*i
) {
1155 i-
>data.val
= i-
>data.val
&-65536;
1158 @ |fraction_to_scaled| rounds a |fraction| and converts it to |scaled|
1160 void mp_fraction_to_round_scaled
(mp_number
*x_orig
) {
1161 int x
= x_orig-
>data.val
;
1162 x_orig-
>type
= mp_scaled_type
;
1163 x_orig-
>data.val
= (x
>=2048 ?
1+((x-2048
) / 4096) : ( x
>=-2048 ?
0 : -(1+((-(x
+1)-2048) / 4096))));
1168 @
* Algebraic and transcendental functions.
1169 \MP\ computes all of the necessary special functions from scratch
, without
1170 relying on |real| arithmetic or system subroutines for sines
, cosines
, etc.
1172 @ To get the square root of a |scaled| number |x|
, we want to calculate
1173 $s
=\lfloor
2^
8\
!\sqrt x
+{1\over2
}\rfloor$. If $x
>0$
, this is the unique
1174 integer such that $
2^
{16}x-s\L s^
2<2^
{16}x
+s$. The following subroutine
1175 determines $s$ by an iterative method that maintains the invariant
1176 relations $x
=2^
{46-2k
}x_0\bmod
2^
{30}$
, $
0<y
=\lfloor
2^
{16-2k
}x_0\rfloor
1177 -s^
2+s\L q
=2s$
, where $x_0$ is the initial value of $x$. The value of~$y$
1178 might
, however
, be zero at the start of the first iteration.
1181 void mp_square_rt
(MP mp
, mp_number
*ret
, mp_number x_orig
) { /* return
, x
: scaled
*/
1183 quarterword k
; /* iteration control counter
*/
1184 integer y
; /* register for intermediate calculations
*/
1185 integer q
; /* register for intermediate calculations
*/
1186 x
= x_orig.data.val
;
1188 @
<Handle square root of zero or negative argument@
>;
1192 while
(x
< fraction_two
) { /* i.e.
, |while x
<@t$
2^
{29}$@
>|\unskip
*/
1196 if
(x
< fraction_four
)
1199 x
= x
- fraction_four
;
1203 @
<Decrease |k| by
1, maintaining the invariant
1204 relations between |x|
, |y|
, and~|q|@
>;
1206 ret-
>data.val
= (int
) (halfp
(q
));
1211 @ @
<Handle square root of zero...@
>=
1215 const char
*hlp
[] = {
1216 "Since I don't take square roots of negative numbers,",
1217 "I'm zeroing this one. Proceed, with fingers crossed.",
1219 mp_snprintf
(msg
, 256, "Square root of %s has been replaced by 0", mp_string_scaled
(mp
, x
));
1220 @.Square root...replaced by
0@
>;
1221 mp_error
(mp
, msg
, hlp
, true
);
1228 @ @
<Decrease |k| by
1, maintaining...@
>=
1231 if
(x
>= fraction_four
) { /* note that |fraction_four
=@t$
2^
{30}$@
>|
*/
1232 x
= x
- fraction_four
;
1238 if
(x
>= fraction_four
) {
1239 x
= x
- fraction_four
;
1245 } else if
(y
<= 0) {
1251 @ Pythagorean addition $\psqrt
{a^
2+b^
2}$ is implemented by an elegant
1252 iterative scheme due to Cleve Moler and Donald Morrison
[{\sl IBM Journal
1253 @^Moler
, Cleve Barry@
>
1254 @^Morrison
, Donald Ross@
>
1255 of Research and Development\
/ \bf27
} (1983), 577--581]. It modifies |a| and~|b|
1256 in such a way that their Pythagorean sum remains invariant
, while the
1257 smaller argument decreases.
1260 void mp_pyth_add
(MP mp
, mp_number
*ret
, mp_number a_orig
, mp_number b_orig
) {
1261 int a
, b
; /* a
,b
: scaled
*/
1262 int r
; /* register used to transform |a| and |b|
, fraction
*/
1263 boolean big
; /* is the result dangerously near $
2^
{31}$?
*/
1264 a
= abs
(a_orig.data.val
);
1265 b
= abs
(b_orig.data.val
);
1270 }; /* now |
0<=b
<=a|
*/
1272 if
(a
< fraction_two
) {
1278 }; /* we reduced the precision to avoid arithmetic overflow
*/
1279 @
<Replace |a| by an approximation to $\psqrt
{a^
2+b^
2}$@
>;
1281 if
(a
< fraction_two
) {
1284 mp-
>arith_error
= true
;
1293 @ The key idea here is to reflect the vector $
(a
,b
)$ about the
1294 line through $
(a
,b
/2)$.
1296 @
<Replace |a| by an approximation to $\psqrt
{a^
2+b^
2}$@
>=
1298 r
= mp_make_fraction
(mp
, b
, a
);
1299 r
= mp_take_fraction
(mp
, r
, r
); /* now $r\approx b^
2/a^
2$
*/
1302 r
= mp_make_fraction
(mp
, r
, fraction_four
+ r
);
1303 a
= a
+ mp_take_fraction
(mp
, a
+ a
, r
);
1304 b
= mp_take_fraction
(mp
, b
, r
);
1308 @ Here is a similar algorithm for $\psqrt
{a^
2-b^
2}$.
1309 It converges slowly when $b$ is near $a$
, but otherwise it works fine.
1312 void mp_pyth_sub
(MP mp
, mp_number
*ret
, mp_number a_orig
, mp_number b_orig
) {
1313 int a
, b
; /* a
,b
: scaled
*/
1314 int r
; /* register used to transform |a| and |b|
, fraction
*/
1315 boolean big
; /* is the result dangerously near $
2^
{31}$?
*/
1316 a
= abs
(a_orig.data.val
);
1317 b
= abs
(b_orig.data.val
);
1319 @
<Handle erroneous |pyth_sub| and set |a
:=0|@
>;
1321 if
(a
< fraction_four
) {
1324 a
= (integer
) halfp
(a
);
1325 b
= (integer
) halfp
(b
);
1328 @
<Replace |a| by an approximation to $\psqrt
{a^
2-b^
2}$@
>;
1336 @ @
<Replace |a| by an approximation to $\psqrt
{a^
2-b^
2}$@
>=
1338 r
= mp_make_fraction
(mp
, b
, a
);
1339 r
= mp_take_fraction
(mp
, r
, r
); /* now $r\approx b^
2/a^
2$
*/
1342 r
= mp_make_fraction
(mp
, r
, fraction_four
- r
);
1343 a
= a
- mp_take_fraction
(mp
, a
+ a
, r
);
1344 b
= mp_take_fraction
(mp
, b
, r
);
1348 @ @
<Handle erroneous |pyth_sub| and set |a
:=0|@
>=
1352 const char
*hlp
[] = {
1353 "Since I don't take square roots of negative numbers,",
1354 "I'm zeroing this one. Proceed, with fingers crossed.",
1356 char
*astr
= strdup
(mp_string_scaled
(mp
, a
));
1358 mp_snprintf
(msg
, 256, "Pythagorean subtraction %s+-+%s has been replaced by 0", astr
, mp_string_scaled
(mp
, b
));
1361 mp_error
(mp
, msg
, hlp
, true
);
1367 @ The subroutines for logarithm and exponential involve two tables.
1368 The first is simple
: |two_to_the
[k
]| equals $
2^k$. The second involves
1369 a bit more calculation
, which the author claims to have done correctly
:
1370 |spec_log
[k
]| is $
2^
{27}$ times $\ln\bigl
(1/(1-2^
{-k
})\bigr
)=
1371 2^
{-k
}+{1\over2
}2^
{-2k
}+{1\over3
}2^
{-3k
}+\cdots\
,$
, rounded to the
1374 @d two_to_the
(A
) (1<<(unsigned
)(A
))
1377 static const integer spec_log
[29] = { 0, /* special logarithms
*/
1378 93032640, 38612034, 17922280, 8662214, 4261238, 2113709,
1379 1052693, 525315, 262400, 131136, 65552, 32772, 16385,
1380 8192, 4096, 2048, 1024, 512, 256, 128, 64, 32, 16, 8, 4, 2, 1, 1
1384 @ Here is the routine that calculates $
2^
8$ times the natural logarithm
1385 of a |scaled| quantity
; it is an integer approximation to $
2^
{24}\ln
(x
/2^
{16})$
,
1386 when |x| is a given positive integer.
1388 The method is based on exercise
1.2.2--25 in
{\sl The Art of Computer
1389 Programming\
/}: During the main iteration we have $
1\L
2^
{-30}x
<1/(1-2^
{1-k
})$
,
1390 and the logarithm of $
2^
{30}x$ remains to be added to an accumulator
1391 register called~$y$. Three auxiliary bits of accuracy are retained in~$y$
1392 during the calculation
, and sixteen auxiliary bits to extend |y| are
1393 kept in~|z| during the initial argument reduction.
(We add
1394 $
100\cdot2^
{16}=6553600$ to~|z| and subtract
100 from~|y| so that |z| will
1395 not become negative
; also
, the actual amount subtracted from~|y| is~
96,
1396 not
~100, because we want to add~
4 for rounding before the final division by~
8.
)
1399 void mp_m_log
(MP mp
, mp_number
*ret
, mp_number x_orig
) { /* return
, x
: scaled
*/
1401 integer y
, z
; /* auxiliary registers
*/
1402 integer k
; /* iteration counter
*/
1403 x
= x_orig.data.val
;
1405 @
<Handle non-positive logarithm@
>;
1407 y
= 1302456956 + 4 - 100; /* $
14\times2^
{27}\ln2\approx1302456956.421063$
*/
1408 z
= 27595 + 6553600; /* and $
2^
{16}\times
.421063\approx
27595$
*/
1409 while
(x
< fraction_four
) {
1413 } /* $
2^
{27}\ln2\approx
93032639.74436163$ and $
2^
{16}\times
.74436163\approx
48782$
*/
1414 y
= y
+ (z
/ unity
);
1416 while
(x
> fraction_four
+ 4) {
1417 @
<Increase |k| until |x| can be multiplied by a
1418 factor of $
2^
{-k
}$
, and adjust $y$ accordingly@
>;
1420 ret-
>data.val
= (y
/ 8);
1425 @ @
<Increase |k| until |x| can...@
>=
1427 z
= ((x
- 1) / two_to_the
(k
)) + 1; /* $z
=\lceil x
/2^k\rceil$
*/
1428 while
(x
< fraction_four
+ z
) {
1437 @ @
<Handle non-positive logarithm@
>=
1440 const char
*hlp
[] = {
1441 "Since I don't take logs of non-positive numbers,",
1442 "I'm zeroing this one. Proceed, with fingers crossed.",
1444 mp_snprintf
(msg
, 256, "Logarithm of %s has been replaced by 0", mp_string_scaled
(mp
, x
));
1445 @.Logarithm...replaced by
0@
>;
1446 mp_error
(mp
, msg
, hlp
, true
);
1451 @ Conversely
, the exponential routine calculates $\exp
(x
/2^
8)$
,
1452 when |x| is |scaled|. The result is an integer approximation to
1453 $
2^
{16}\exp
(x
/2^
{24})$
, when |x| is regarded as an integer.
1456 void mp_m_exp
(MP mp
, mp_number
*ret
, mp_number x_orig
) {
1457 quarterword k
; /* loop control index
*/
1458 integer y
, z
; /* auxiliary registers
*/
1460 x
= x_orig.data.val
;
1461 if
(x
> 174436200) {
1462 /* $
2^
{24}\ln
((2^
{31}-1)/2^
{16})\approx
174436199.51$
*/
1463 mp-
>arith_error
= true
;
1464 ret-
>data.val
= EL_GORDO
;
1465 } else if
(x
< -197694359) {
1466 /* $
2^
{24}\ln
(2^
{-1}/2^
{16})\approx-197694359.45$
*/
1471 y
= 04000000; /* $y
=2^
{20}$
*/
1473 if
(x
<= 127919879) {
1474 z
= 1023359037 - 8 * x
;
1475 /* $
2^
{27}\ln
((2^
{31}-1)/2^
{20})\approx
1023359037.125$
*/
1477 z
= 8 * (174436200 - x
); /* |z| is always nonnegative
*/
1481 @
<Multiply |y| by $\exp
(-z
/2^
{27})$@
>;
1483 ret-
>data.val
= ((y
+ 8) / 16);
1490 @ The idea here is that subtracting |spec_log
[k
]| from |z| corresponds
1491 to multiplying |y| by $
1-2^
{-k
}$.
1493 A subtle point
(which had to be checked
) was that if $x
=127919879$
, the
1494 value of~|y| will decrease so that |y
+8| doesn't overflow. In fact
,
1495 $z$ will be
5 in this case
, and |y| will decrease by~
64 when |k
=25|
1496 and by~
16 when |k
=27|.
1498 @
<Multiply |y| by...@
>=
1501 while
(z
>= spec_log
[k
]) {
1503 y
= y
- 1 - ((y
- two_to_the
(k
- 1)) / two_to_the
(k
));
1508 @ The trigonometric subroutines use an auxiliary table such that
1509 |spec_atan
[k
]| contains an approximation to the |angle| whose tangent
1510 is~$
1/2^k$. $\arctan2^
{-k
}$ times $
2^
{20}\cdot180
/\pi$
1513 static const int spec_atan
[27] = { 0, 27855475, 14718068, 7471121, 3750058,
1514 1876857, 938658, 469357, 234682, 117342, 58671, 29335, 14668, 7334, 3667,
1515 1833, 917, 458, 229, 115, 57, 29, 14, 7, 4, 2, 1
1519 @ Given integers |x| and |y|
, not both zero
, the |n_arg| function
1520 returns the |angle| whose tangent points in the direction $
(x
,y
)$.
1521 This subroutine first determines the correct octant
, then solves the
1522 problem for |
0<=y
<=x|
, then converts the result appropriately to
1523 return an answer in the range |
-one_eighty_deg
<=@t$\theta$@
><=one_eighty_deg|.
1524 (The answer is |
+one_eighty_deg| if |y
=0| and |x
<0|
, but an answer of
1525 |
-one_eighty_deg| is possible if
, for example
, |y
=-1| and $x
=-2^
{30}$.
)
1527 The octants are represented in a ``Gray code
,'' since that turns out
1528 to be computationally simplest.
1534 @d second_octant
(first_octant
+switch_x_and_y
)
1535 @d third_octant
(first_octant
+switch_x_and_y
+negate_x
)
1536 @d fourth_octant
(first_octant
+negate_x
)
1537 @d fifth_octant
(first_octant
+negate_x
+negate_y
)
1538 @d sixth_octant
(first_octant
+switch_x_and_y
+negate_x
+negate_y
)
1539 @d seventh_octant
(first_octant
+switch_x_and_y
+negate_y
)
1540 @d eighth_octant
(first_octant
+negate_y
)
1543 void mp_n_arg
(MP mp
, mp_number
*ret
, mp_number x_orig
, mp_number y_orig
) {
1544 integer z
; /* auxiliary register
*/
1545 integer t
; /* temporary storage
*/
1546 quarterword k
; /* loop counter
*/
1547 int octant
; /* octant code
*/
1549 x
= x_orig.data.val
;
1550 y
= y_orig.data.val
;
1552 octant
= first_octant
;
1555 octant
= first_octant
+ negate_x
;
1559 octant
= octant
+ negate_y
;
1565 octant
= octant
+ switch_x_and_y
;
1568 @
<Handle undefined arg@
>;
1570 ret-
>type
= mp_angle_type
;
1571 @
<Set variable |z| to the arg of $
(x
,y
)$@
>;
1572 @
<Return an appropriate answer based on |z| and |octant|@
>;
1577 @ @
<Handle undefined arg@
>=
1579 const char
*hlp
[] = {
1580 "The `angle' between two identical points is undefined.",
1581 "I'm zeroing this one. Proceed, with fingers crossed.",
1583 mp_error
(mp
, "angle(0,0) is taken as zero", hlp
, true
);
1584 @.angle
(0,0)...zero@
>;
1589 @ @
<Return an appropriate answer...@
>=
1595 ret-
>data.val
= (ninety_deg
- z
);
1598 ret-
>data.val
= (ninety_deg
+ z
);
1601 ret-
>data.val
= (one_eighty_deg
- z
);
1604 ret-
>data.val
= (z
- one_eighty_deg
);
1607 ret-
>data.val
= (-z
- ninety_deg
);
1609 case seventh_octant
:
1610 ret-
>data.val
= (z
- ninety_deg
);
1613 ret-
>data.val
= (-z
);
1615 } /* there are no other cases
*/
1618 @ At this point we have |x
>=y
>=0|
, and |x
>0|. The numbers are scaled up
1619 or down until $
2^
{28}\L x
<2^
{29}$
, so that accurate fixed-point calculations
1622 @
<Set variable |z| to the arg...@
>=
1623 while
(x
>= fraction_two
) {
1629 while
(x
< fraction_one
) {
1633 @
<Increase |z| to the arg of $
(x
,y
)$@
>;
1636 @ During the calculations of this section
, variables |x| and~|y|
1637 represent actual coordinates $
(x
,2^
{-k
}y
)$. We will maintain the
1638 condition |x
>=y|
, so that the tangent will be at most $
2^
{-k
}$.
1639 If $x
<2y$
, the tangent is greater than $
2^
{-k-1
}$. The transformation
1640 $
(a
,b
)\mapsto
(a
+b\tan\phi
,b-a\tan\phi
)$ replaces $
(a
,b
)$ by
1641 coordinates whose angle has decreased by~$\phi$
; in the special case
1642 $a
=x$
, $b
=2^
{-k
}y$
, and $\tan\phi
=2^
{-k-1
}$
, this operation reduces
1643 to the particularly simple iteration shown here.
[Cf.~John E. Meggitt
,
1644 @^Meggitt
, John E.@
>
1645 {\sl IBM Journal of Research and Development\
/ \bf6
} (1962), 210--226.
]
1647 The initial value of |x| will be multiplied by at most
1648 $
(1+{1\over2
})(1+{1\over8
})(1+{1\over32
})\cdots\approx
1.7584$
; hence
1649 there is no chance of integer overflow.
1651 @
<Increase |z|...@
>=
1657 z
= z
+ spec_atan
[k
];
1659 x
= x
+ (y
/ two_to_the
(k
+ k
));
1667 z
= z
+ spec_atan
[k
];
1672 @ Conversely
, the |n_sin_cos| routine takes an |angle| and produces the sine
1673 and cosine of that angle. The results of this routine are
1674 stored in global integer variables |n_sin| and |n_cos|.
1676 @ Given an integer |z| that is $
2^
{20}$ times an angle $\theta$ in degrees
,
1677 the purpose of |n_sin_cos
(z
)| is to set
1678 |x
=@t$r\cos\theta$@
>| and |y
=@t$r\sin\theta$@
>|
(approximately
),
1679 for some rather large number~|r|. The maximum of |x| and |y|
1680 will be between $
2^
{28}$ and $
2^
{30}$
, so that there will be hardly
1681 any loss of accuracy. Then |x| and~|y| are divided by~|r|.
1683 @d forty_five_deg
0264000000 /* $
45\cdot2^
{20}$
, represents $
45^\circ$
*/
1684 @d ninety_deg
0550000000 /* $
90\cdot2^
{20}$
, represents $
90^\circ$
*/
1685 @d one_eighty_deg
01320000000 /* $
180\cdot2^
{20}$
, represents $
180^\circ$
*/
1686 @d three_sixty_deg
02640000000 /* $
360\cdot2^
{20}$
, represents $
360^\circ$
*/
1688 @d odd
(A
) (abs
(A
)%2==1)
1690 @ Compute a multiple of the sine and cosine
1693 void mp_n_sin_cos
(MP mp
, mp_number z_orig
, mp_number
*n_cos
, mp_number
*n_sin
) {
1694 quarterword k
; /* loop control variable
*/
1695 int q
; /* specifies the quadrant
*/
1696 integer x
, y
, t
; /* temporary registers
*/
1698 mp_number x_n
, y_n
, ret
;
1702 z
= z_orig.data.val
;
1704 z
= z
+ three_sixty_deg
;
1705 z
= z
% three_sixty_deg
; /* now |
0<=z
<three_sixty_deg|
*/
1706 q
= z
/ forty_five_deg
;
1707 z
= z
% forty_five_deg
;
1711 z
= forty_five_deg
- z
;
1712 @
<Subtract angle |z| from |
(x
,y
)|@
>;
1713 @
<Convert |
(x
,y
)| to the octant determined by~|q|@
>;
1716 mp_pyth_add
(mp
, &ret, x_n, y_n);
1717 n_cos-
>data.val
= mp_make_fraction
(mp
, x
, ret.data.val
);
1718 n_sin-
>data.val
= mp_make_fraction
(mp
, y
, ret.data.val
);
1725 @ In this case the octants are numbered sequentially.
1727 @
<Convert |
(x
,...@
>=
1761 } /* there are no other cases
*/
1764 @ The main iteration of |n_sin_cos| is similar to that of |n_arg| but
1765 applied in reverse. The values of |spec_atan
[k
]| decrease slowly enough
1766 that this loop is guaranteed to terminate before the
(nonexistent
) value
1767 |spec_atan
[27]| would be required.
1769 @
<Subtract angle |z|...@
>=
1772 if
(z
>= spec_atan
[k
]) {
1773 z
= z
- spec_atan
[k
];
1775 x
= t
+ y
/ two_to_the
(k
);
1776 y
= y
- t
/ two_to_the
(k
);
1781 y
= 0 /* this precaution may never be needed
*/
1784 @ To initialize the |randoms| table
, we call the following routine.
1787 void mp_init_randoms
(MP mp
, int seed
) {
1788 int j
, jj
, k
; /* more or less random integers
*/
1789 int i
; /* index into |randoms|
*/
1791 while
(j
>= fraction_one
) {
1795 for
(i
= 0; i
<= 54; i
++) {
1801 mp-
>randoms
[(i
* 21) % 55].data.val
= j
;
1803 mp_new_randoms
(mp
);
1804 mp_new_randoms
(mp
);
1805 mp_new_randoms
(mp
); /* ``warm up'' the array
*/
1810 void mp_print_number
(MP mp
, mp_number n
) {
1811 mp_print_scaled
(mp
, n.data.val
);
1816 char
* mp_number_tostring
(MP mp
, mp_number n
) {
1817 return mp_string_scaled
(mp
, n.data.val
);
1821 void mp_number_modulo
(mp_number
*a
, mp_number b
) {
1822 a-
>data.val
= a-
>data.val
% b.data.val
;
1829 @ To consume a random fraction
, the program below will say `|next_random|'.
1832 static void mp_next_random
(MP mp
, mp_number
*ret
) {
1833 if
( mp-
>j_random
==0 )
1836 mp-
>j_random
= mp-
>j_random-1
;
1837 mp_number_clone
(ret
, mp-
>randoms
[mp-
>j_random
]);
1841 @ To produce a uniform random number in the range |
0<=u
<x| or |
0>=u
>x|
1842 or |
0=u
=x|
, given a |scaled| value~|x|
, we proceed as shown here.
1844 Note that the call of |take_fraction| will produce the values
0 and~|x|
1845 with about half the probability that it will produce any other particular
1846 values between
0 and~|x|
, because it rounds its answers.
1849 static void mp_m_unif_rand
(MP mp
, mp_number
*ret
, mp_number x_orig
) {
1850 mp_number y
; /* trial value
*/
1857 mp_number_clone
(&x, x_orig);
1858 mp_number_clone
(&abs_x, x);
1859 mp_number_abs
(&abs_x);
1860 mp_next_random
(mp
, &u);
1861 /*take_fraction
(y
, abs_x
, u
);*/
1862 mp_number_take_fraction
(mp
,&y, abs_x,u);
1864 if
(mp_number_equal
(y
, abs_x
)) {
1865 /*set_number_to_zero
(*ret
);*/
1866 mp_number_clone
(ret
, ((math_data
*)mp-
>math
)->zero_t
);
1867 } else if
(mp_number_greater
(x
, ((math_data
*)mp-
>math
)->zero_t
)) {
1868 mp_number_clone
(ret
, y
);
1870 mp_number_clone
(ret
, y
);
1871 mp_number_negate
(ret
);
1873 free_number
(abs_x
);
1881 @ Finally
, a normal deviate with mean zero and unit standard deviation
1882 can readily be obtained with the ratio method
(Algorithm
3.4.1R in
1883 {\sl The Art of Computer Programming\
/}).
1886 static void mp_m_norm_rand
(MP mp
, mp_number
*ret
) {
1892 new_number
(ab_vs_cd
);
1902 mp_next_random
(mp
, &v);
1903 mp_number_substract
(&v, ((math_data *)mp->math)->fraction_half_t);
1904 mp_number_take_fraction
(mp
,&xa, ((math_data *)mp->math)->sqrt_8_e_k, v);
1906 mp_next_random
(mp
, &u);
1907 mp_number_clone
(&abs_x, xa);
1908 mp_number_abs
(&abs_x);
1909 } while
(!mp_number_less
(abs_x
, u
));
1910 mp_number_make_fraction
(mp
, &r, xa, u);
1911 mp_number_clone
(&xa, r);
1912 mp_m_log
(mp
,&la, u);
1913 mp_set_number_from_substraction
(&la, ((math_data *)mp->math)->twelve_ln_2_k, la);
1914 mp_ab_vs_cd
(mp
,&ab_vs_cd, ((math_data *)mp->math)->one_k, la, xa, xa);
1915 } while
(mp_number_less
(ab_vs_cd
,((math_data
*)mp-
>math
)->zero_t
));
1916 mp_number_clone
(ret
, xa
);
1917 free_number
(ab_vs_cd
);
1919 free_number
(abs_x
);