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