1
5
6
10
11# include "mpmathposit.h"
12
13# define mp_fraction_multiplier 4096
14# define mp_angle_multiplier 16
15# define mp_warning_limit pow(2.0,52)
16# define odd(A) (abs(A)%2==1)
17# define two_to_the(A) (1<<(unsigned)(A))
18# define set_cur_cmd(A) mp->cur_mod_->command = (A)
19# define set_cur_mod(A) mp->cur_mod_->data.n.data.pval = (A)
20
21typedef struct mp_posit_info {
22 posit_t unity;
23 posit_t zero;
24 posit_t one;
25 posit_t two;
26 posit_t three;
27 posit_t four;
28 posit_t five;
29 posit_t eight;
30 posit_t seven;
31 posit_t sixteen;
32 posit_t half_unit;
33 posit_t minusone;
34 posit_t three_quarter_unit;
35 posit_t d16;
36 posit_t d64;
37 posit_t d256;
38 posit_t d4096;
39 posit_t d65536;
40 posit_t dp90;
41 posit_t dp180;
42 posit_t dp270;
43 posit_t dp360;
44 posit_t dm90;
45 posit_t dm180;
46 posit_t dm270;
47 posit_t dm360;
48 posit_t fraction_multiplier;
49 posit_t negative_fraction_multiplier;
50 posit_t angle_multiplier;
51 posit_t fraction_one;
52 posit_t fraction_two;
53 posit_t fraction_three;
54 posit_t fraction_four;
55 posit_t fraction_half;
56 posit_t fraction_one_and_half;
57 posit_t one_eighty_degrees;
58 posit_t negative_one_eighty_degrees;
59 posit_t three_sixty_degrees;
60 posit_t no_crossing;
61 posit_t one_crossing;
62 posit_t zero_crossing;
63 posit_t error_correction;
64 posit_t pi;
65 posit_t pi_divided_by_180;
66 posit_t epsilon;
67 posit_t EL_GORDO;
68 posit_t negative_EL_GORDO;
69 posit_t one_third_EL_GORDO;
70 posit_t coef;
71 posit_t coef_bound;
72 posit_t scaled_threshold;
73 posit_t fraction_threshold;
74 posit_t equation_threshold;
75 posit_t near_zero_angle;
76 posit_t p_over_v_threshold;
77 posit_t warning_limit;
78 posit_t sqrt_two_mul_fraction_one;
79 posit_t sqrt_five_minus_one_mul_fraction_one_and_half;
80 posit_t three_minus_sqrt_five_mul_fraction_one_and_half;
81 posit_t d180_divided_by_pi_mul_angle;
82 int initialized;
83} mp_posit_info;
84
85static mp_posit_info mp_posit_data = {
86 .initialized = 0,
87};
88
89static inline posit_t mp_posit_aux_make_fraction (posit_t p, posit_t q) { return posit_mul(posit_div(p,q), mp_posit_data.fraction_multiplier); }
90static inline posit_t mp_posit_aux_take_fraction (posit_t p, posit_t q) { return posit_div(posit_mul(p,q), mp_posit_data.fraction_multiplier); }
91static inline posit_t mp_posit_aux_make_scaled (posit_t p, posit_t q) { return posit_div(p,q); }
92
93
97
98static void mp_posit_allocate_number(MP mp, mp_number *n, mp_number_type t)
99{
100 (void) mp;
101 n->data.pval = mp_posit_data.zero;
102 n->type = t;
103}
104
105static void mp_posit_allocate_clone(MP mp, mp_number *n, mp_number_type t, mp_number *v)
106{
107 (void) mp;
108 n->type = t;
109 n->data.pval = v->data.pval;
110}
111
112static void mp_posit_allocate_abs(MP mp, mp_number *n, mp_number_type t, mp_number *v)
113{
114 (void) mp;
115 n->type = t;
116 n->data.pval = posit_fabs(v->data.pval);
117}
118
119static void mp_posit_allocate_double(MP mp, mp_number *n, double v)
120{
121 (void) mp;
122 n->type = mp_scaled_type;
123 n->data.pval = double_to_posit(v);
124}
125
126static void mp_posit_free_number(MP mp, mp_number *n)
127{
128 (void) mp;
129 n->type = mp_nan_type;
130}
131
132static void mp_posit_set_from_int(mp_number *A, int B)
133{
134 A->data.pval = integer_to_posit(B);
135}
136
137static void mp_posit_set_from_boolean(mp_number *A, int B)
138{
139 A->data.pval = integer_to_posit(B);
140}
141
142static void mp_posit_set_from_scaled(mp_number *A, int B)
143{
144 A->data.pval = posit_div(integer_to_posit(B), mp_posit_data.d65536);
145}
146
147static void mp_posit_set_from_double(mp_number *A, double B)
148{
149 A->data.pval = double_to_posit(B);
150}
151
152static void mp_posit_set_from_addition(mp_number *A, mp_number *B, mp_number *C)
153{
154 A->data.pval = posit_add(B->data.pval, C->data.pval);
155}
156
157static void mp_posit_set_half_from_addition(mp_number *A, mp_number *B, mp_number *C)
158{
159 A->data.pval = posit_div(posit_add(B->data.pval,C->data.pval), mp_posit_data.two);
160}
161
162static void mp_posit_set_from_subtraction(mp_number *A, mp_number *B, mp_number *C)
163{
164 A->data.pval = posit_sub(B->data.pval, C->data.pval);
165}
166
167static void mp_posit_set_half_from_subtraction(mp_number *A, mp_number *B, mp_number *C)
168{
169 A->data.pval = posit_div(posit_sub(B->data.pval, C->data.pval), mp_posit_data.two);
170}
171
172static void mp_posit_set_from_div(mp_number *A, mp_number *B, mp_number *C)
173{
174 A->data.pval = posit_div(B->data.pval, C->data.pval);
175}
176
177static void mp_posit_set_from_mul(mp_number *A, mp_number *B, mp_number *C)
178{
179 A->data.pval = posit_mul(B->data.pval, C->data.pval);
180}
181
182static void mp_posit_set_from_int_div(mp_number *A, mp_number *B, int C)
183{
184 A->data.pval = posit_div(B->data.pval, integer_to_posit(C));
185}
186
187static void mp_posit_set_from_int_mul(mp_number *A, mp_number *B, int C)
188{
189 A->data.pval = posit_mul(B->data.pval, integer_to_posit(C));
190}
191
192static void mp_posit_set_from_of_the_way(MP mp, mp_number *A, mp_number *t, mp_number *B, mp_number *C)
193{
194 (void) mp;
195 A->data.pval = posit_sub(B->data.pval, mp_posit_aux_take_fraction(posit_sub(B->data.pval, C->data.pval), t->data.pval));
196}
197
198static void mp_posit_negate(mp_number *A)
199{
200 A->data.pval = posit_neg(A->data.pval);
201}
202
203static void mp_posit_add(mp_number *A, mp_number *B)
204{
205 A->data.pval = posit_add(A->data.pval, B->data.pval);
206}
207
208static void mp_posit_subtract(mp_number *A, mp_number *B)
209{
210 A->data.pval = posit_sub(A->data.pval, B->data.pval);
211}
212
213static void mp_posit_half(mp_number *A)
214{
215 A->data.pval = posit_div(A->data.pval, mp_posit_data.two);
216}
217
218static void mp_posit_double(mp_number *A)
219{
220 A->data.pval = posit_mul(A->data.pval, mp_posit_data.two);
221}
222
223static void mp_posit_add_scaled(mp_number *A, int B)
224{
225
226 A->data.pval = posit_add(A->data.pval, posit_div(integer_to_posit(B), mp_posit_data.d65536));
227}
228
229static void mp_posit_multiply_int(mp_number *A, int B)
230{
231 A->data.pval = posit_mul(A->data.pval, integer_to_posit(B));
232}
233
234static void mp_posit_divide_int(mp_number *A, int B)
235{
236 A->data.pval = posit_div(A->data.pval, integer_to_posit(B));
237}
238
239static void mp_posit_abs(mp_number *A)
240{
241 A->data.pval = posit_fabs(A->data.pval);
242}
243
244static void mp_posit_clone(mp_number *A, mp_number *B)
245{
246 A->data.pval = B->data.pval;
247}
248
249static void mp_posit_negated_clone(mp_number *A, mp_number *B)
250{
251 A->data.pval = posit_neg(B->data.pval);
252}
253
254static void mp_posit_abs_clone(mp_number *A, mp_number *B)
255{
256 A->data.pval = posit_fabs(B->data.pval);
257}
258
259static void mp_posit_swap(mp_number *A, mp_number *B)
260{
261 posit_t swap_tmp = A->data.pval;
262 A->data.pval = B->data.pval;
263 B->data.pval = swap_tmp;
264}
265
266static void mp_posit_fraction_to_scaled(mp_number *A)
267{
268 A->type = mp_scaled_type;
269 A->data.pval = posit_div(A->data.pval, mp_posit_data.fraction_multiplier);
270}
271
272static void mp_posit_angle_to_scaled(mp_number *A)
273{
274 A->type = mp_scaled_type;
275 A->data.pval = posit_div(A->data.pval, mp_posit_data.angle_multiplier);
276}
277
278static void mp_posit_scaled_to_fraction(mp_number *A)
279{
280 A->type = mp_fraction_type;
281 A->data.pval = posit_mul(A->data.pval, mp_posit_data.fraction_multiplier);
282}
283
284static void mp_posit_scaled_to_angle(mp_number *A)
285{
286 A->type = mp_angle_type;
287 A->data.pval = posit_mul(A->data.pval, mp_posit_data.angle_multiplier);
288}
289
290static int mp_posit_to_scaled(mp_number *A)
291{
292 return (int) posit_to_integer(posit_mul(A->data.pval, mp_posit_data.d65536));
293}
294
295static int mp_posit_to_int(mp_number *A)
296{
297 return (int) posit_to_integer(A->data.pval);
298}
299
300static int mp_posit_to_boolean(mp_number *A)
301{
302 return posit_eq_zero(A->data.pval) ? 0 : 1;
303}
304
305static double mp_posit_to_double(mp_number *A)
306{
307 return posit_to_double(A->data.pval);
308}
309
310static int mp_posit_odd(mp_number *A)
311{
312 return (int) odd(posit_to_integer(A->data.pval));
313}
314
315static int mp_posit_equal(mp_number *A, mp_number *B)
316{
317 return posit_eq(A->data.pval, B->data.pval);
318}
319
320static int mp_posit_greater(mp_number *A, mp_number *B)
321{
322 return posit_gt(A->data.pval, B->data.pval);
323}
324
325static int mp_posit_less(mp_number *A, mp_number *B)
326{
327 return posit_lt(A->data.pval, B->data.pval);
328}
329
330static int mp_posit_non_equal_abs(mp_number *A, mp_number *B)
331{
332 return ! posit_eq(posit_fabs(A->data.pval), posit_fabs(B->data.pval));
333}
334
335static char *mp_posit_number_tostring(MP mp, mp_number *n)
336{
337 static char set[64];
338 int l = 0;
339 char *ret = mp_memory_allocate(64);
340 (void) mp;
341 snprintf(set, 64, "%.20g", posit_to_double(n->data.pval));
342 while (set[l] == ' ') {
343 l++;
344 }
345 strcpy(ret, set+l);
346 return ret;
347}
348
349static void mp_posit_print_number(MP mp, mp_number *n)
350{
351 char *str = mp_posit_number_tostring(mp, n);
352 mp_print_e_str(mp, str);
353 mp_memory_free(str);
354}
355
356
357
358static void mp_posit_slow_add(MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_orig)
359{
360 if (posit_gt(x_orig->data.pval, mp_posit_data.zero)) {
361 if (posit_le(y_orig->data.pval, posit_sub(mp_posit_data.EL_GORDO, x_orig->data.pval))) {
362 ret->data.pval = posit_add(x_orig->data.pval, y_orig->data.pval);
363 } else {
364 mp->arith_error = 1;
365 ret->data.pval = mp_posit_data.EL_GORDO;
366 }
367 } else if (posit_le(posit_neg(y_orig->data.pval), posit_add(mp_posit_data.EL_GORDO, x_orig->data.pval))) {
368 ret->data.pval = posit_add(x_orig->data.pval, y_orig->data.pval);
369 } else {
370 mp->arith_error = 1;
371 ret->data.pval = mp_posit_data.negative_EL_GORDO;
372 }
373}
374
375static void mp_posit_make_fraction(MP mp, mp_number *ret, mp_number *p, mp_number *q) {
376 (void) mp;
377 ret->data.pval = mp_posit_aux_make_fraction(p->data.pval, q->data.pval);
378}
379
380static void mp_posit_take_fraction(MP mp, mp_number *ret, mp_number *p, mp_number *q) {
381 (void) mp;
382 ret->data.pval = mp_posit_aux_take_fraction(p->data.pval, q->data.pval);
383}
384
385static void mp_posit_take_scaled(MP mp, mp_number *ret, mp_number *p_orig, mp_number *q_orig)
386{
387 (void) mp;
388 ret->data.pval = posit_mul(p_orig->data.pval, q_orig->data.pval);
389}
390
391static void mp_posit_make_scaled(MP mp, mp_number *ret, mp_number *p_orig, mp_number *q_orig)
392{
393 (void) mp;
394 ret->data.pval = posit_div(p_orig->data.pval, q_orig->data.pval);
395}
396
397static void mp_posit_aux_wrapup_numeric_token(MP mp, unsigned char *start, unsigned char *stop)
398{
399 double result;
400 char *end = (char *) stop;
401 errno = 0;
402 result = strtod((char *) start, &end);
403 if (errno == 0) {
404 set_cur_mod(double_to_posit(result));
405 if (result >= mp_warning_limit) {
406 if (posit_gt(internal_value(mp_warning_check_internal).data.pval, mp_posit_data.zero) && (mp->scanner_status != mp_tex_flushing_state)) {
407 char msg[256];
408 snprintf(msg, 256, "Number is too large (%g)", result);
409 mp_error(
410 mp,
411 msg,
412 "Continue and I'll try to cope with that big value; but it might be dangerous."
413 "(Set warningcheck := 0 to suppress this message.)"
414 );
415 }
416 }
417 } else if (mp->scanner_status != mp_tex_flushing_state) {
418 mp_error(
419 mp,
420 "Enormous number has been reduced.",
421 "I could not handle this number specification probably because it is out of"
422 "range."
423 );
424 set_cur_mod(mp_posit_data.EL_GORDO);
425 }
426 set_cur_cmd(mp_numeric_command);
427}
428
429static void mp_posit_aux_find_exponent(MP mp)
430{
431 if (mp->buffer[mp->cur_input.loc_field] == 'e' || mp->buffer[mp->cur_input.loc_field] == 'E') {
432 mp->cur_input.loc_field++;
433 if (!(mp->buffer[mp->cur_input.loc_field] == '+'
434 || mp->buffer[mp->cur_input.loc_field] == '-'
435 || mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class)) {
436 mp->cur_input.loc_field--;
437 return;
438 }
439 if (mp->buffer[mp->cur_input.loc_field] == '+'
440 || mp->buffer[mp->cur_input.loc_field] == '-') {
441 mp->cur_input.loc_field++;
442 }
443 while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class) {
444 mp->cur_input.loc_field++;
445 }
446 }
447}
448
449static void mp_posit_scan_fractional_token(MP mp, int n)
450{
451 unsigned char *start = &mp->buffer[mp->cur_input.loc_field -1];
452 unsigned char *stop;
453 (void) n;
454 while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class) {
455 mp->cur_input.loc_field++;
456 }
457 mp_posit_aux_find_exponent(mp);
458 stop = &mp->buffer[mp->cur_input.loc_field-1];
459 mp_posit_aux_wrapup_numeric_token(mp, start, stop);
460}
461
462static void mp_posit_scan_numeric_token(MP mp, int n)
463{
464 unsigned char *start = &mp->buffer[mp->cur_input.loc_field -1];
465 unsigned char *stop;
466 (void) n;
467 while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class) {
468 mp->cur_input.loc_field++;
469 }
470 if (mp->buffer[mp->cur_input.loc_field] == '.' && mp->buffer[mp->cur_input.loc_field+1] != '.') {
471 mp->cur_input.loc_field++;
472 while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class) {
473 mp->cur_input.loc_field++;
474 }
475 }
476 mp_posit_aux_find_exponent(mp);
477 stop = &mp->buffer[mp->cur_input.loc_field-1];
478 mp_posit_aux_wrapup_numeric_token(mp, start, stop);
479}
480
481static void mp_posit_velocity(MP mp, mp_number *ret, mp_number *st, mp_number *ct, mp_number *sf, mp_number *cf, mp_number *t)
482{
483 posit_t acc, num, denom;
484 (void) mp;
485 acc = mp_posit_aux_take_fraction(
486 mp_posit_aux_take_fraction(
487 posit_sub(st->data.pval, posit_div(sf->data.pval, mp_posit_data.sixteen)),
488 posit_sub(sf->data.pval, posit_div(st->data.pval, mp_posit_data.sixteen))
489 ),
490 posit_sub(ct->data.pval,cf->data.pval)
491 );
492 num = posit_add(
493 mp_posit_data.fraction_two,
494 mp_posit_aux_take_fraction(
495 acc,
496 mp_posit_data.sqrt_two_mul_fraction_one
497 )
498 );
499 denom = posit_add(
500 mp_posit_data.fraction_three,
501 posit_add(
502 mp_posit_aux_take_fraction(
503 ct->data.pval,
504 mp_posit_data.sqrt_five_minus_one_mul_fraction_one_and_half
505 ),
506 mp_posit_aux_take_fraction(
507 cf->data.pval,
508 mp_posit_data.three_minus_sqrt_five_mul_fraction_one_and_half
509 )
510 )
511 );
512 if (posit_ne(t->data.pval, mp_posit_data.unity)) {
513 num = mp_posit_aux_make_scaled(num, t->data.pval);
514 }
515 if (posit_ge(posit_div(num, mp_posit_data.four), denom)) {
516 ret->data.pval = mp_posit_data.fraction_four;
517 } else {
518 ret->data.pval = mp_posit_aux_make_fraction(num, denom);
519 }
520}
521
522static int mp_posit_ab_vs_cd (mp_number *a_orig, mp_number *b_orig, mp_number *c_orig, mp_number *d_orig)
523{
524 posit_t ab = posit_mul(a_orig->data.pval, b_orig->data.pval);
525 posit_t cd = posit_mul(c_orig->data.pval, d_orig->data.pval);
526 if (posit_eq(ab,cd)) {
527 return 0;
528 } else if (posit_lt(ab,cd)) {
529 return -1;
530 } else {
531 return 1;
532 }
533}
534
535static void mp_posit_crossing_point(MP mp, mp_number *ret, mp_number *aa, mp_number *bb, mp_number *cc)
536{
537 posit_t d;
538 posit_t xx, x0, x1, x2;
539 posit_t a = aa->data.pval;
540 posit_t b = bb->data.pval;
541 posit_t c = cc->data.pval;
542 (void) mp;
543 if (posit_lt(a, mp_posit_data.zero)) {
544 ret->data.pval = mp_posit_data.zero_crossing;
545 return;
546 }
547 if (posit_ge(c, mp_posit_data.zero)) {
548 if (posit_ge(b, mp_posit_data.zero)) {
549 if (posit_gt(c, mp_posit_data.zero)) {
550 ret->data.pval = mp_posit_data.no_crossing;
551 } else if (posit_eq_zero(a) && posit_eq_zero(b)) {
552 ret->data.pval = mp_posit_data.no_crossing;
553 } else {
554 ret->data.pval = mp_posit_data.one_crossing;
555 }
556 return;
557 }
558 if (posit_eq_zero(a)) {
559 ret->data.pval = mp_posit_data.zero_crossing;
560 return;
561 }
562 } else if (posit_eq_zero(a) && posit_le(b, mp_posit_data.zero)) {
563 ret->data.pval = mp_posit_data.zero_crossing;
564 return;
565 }
566
567 d = mp_posit_data.epsilon;
568 x0 = a;
569 x1 = posit_sub(a, b);
570 x2 = posit_sub(b, c);
571 do {
572
573 posit_t x = posit_add(posit_div(posit_add(x1, x2), mp_posit_data.two), mp_posit_data.error_correction);
574 if (posit_gt(posit_sub(x1, x0), x0)) {
575 x2 = x;
576 x0 = posit_add(x0, x0);
577 d = posit_add(d, d);
578 } else {
579 xx = posit_sub(posit_add(x1, x), x0);
580 if (posit_gt(xx, x0)) {
581 x2 = x;
582 x0 = posit_add(x0, x0);
583 d = posit_add(d, d);
584 } else {
585 x0 = posit_sub(x0, xx);
586 if (posit_le(x, x0) && posit_le(posit_add(x, x2), x0)) {
587 ret->data.pval = mp_posit_data.no_crossing;
588 return;
589 }
590 x1 = x;
591 d = posit_add(posit_add(d, d), mp_posit_data.epsilon);
592 }
593 }
594 } while (posit_lt(d, mp_posit_data.fraction_one));
595 ret->data.pval = posit_sub(d, mp_posit_data.fraction_one);
596}
597
598
599
600static int mp_posit_round_unscaled(mp_number *x_orig)
601{
602 return posit_i_round(x_orig->data.pval);
603}
604
605static void mp_posit_floor(mp_number *i)
606{
607 i->data.pval = posit_floor(i->data.pval);
608}
609
610static void mp_posit_fraction_to_round_scaled(mp_number *x_orig)
611{
612 x_orig->type = mp_scaled_type;
613 x_orig->data.pval = posit_div(x_orig->data.pval, mp_posit_data.fraction_multiplier);
614}
615
616static void mp_posit_square_rt(MP mp, mp_number *ret, mp_number *x_orig)
617{
618 if (posit_gt(x_orig->data.pval, mp_posit_data.zero)) {
619 ret->data.pval = posit_sqrt(x_orig->data.pval);
620 } else {
621 if (posit_lt(x_orig->data.pval, mp_posit_data.zero)) {
622 char msg[256];
623 char *xstr = mp_posit_number_tostring(mp, x_orig);
624 snprintf(msg, 256, "Square root of %s has been replaced by 0", xstr);
625 mp_memory_free(xstr);
626 mp_error(
627 mp,
628 msg,
629 "Since I don't take square roots of negative numbers, I'm zeroing this one.\n"
630 "Proceed, with fingers crossed."
631 );
632 }
633 ret->data.pval = mp_posit_data.zero;
634 }
635}
636
637static void mp_posit_pyth_add(MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig)
638{
639 (void) mp;
640 ret->data.pval = posit_sqrt(
641 posit_add(
642 posit_mul(
643 a_orig->data.pval,
644 a_orig->data.pval
645 ),
646 posit_mul(
647 b_orig->data.pval,
648 b_orig->data.pval
649 )
650 )
651 );
652}
653
654static void mp_posit_pyth_sub(MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig)
655{
656
657 if (posit_gt(a_orig->data.pval,b_orig->data.pval)) {
658 a_orig->data.pval = posit_sqrt(
659 posit_sub(
660 posit_mul(
661 a_orig->data.pval,
662 a_orig->data.pval
663 ),
664 posit_mul(
665 b_orig->data.pval,
666 b_orig->data.pval
667 )
668 )
669 );
670 } else {
671 if (posit_lt(a_orig->data.pval,b_orig->data.pval)) {
672 char msg[256];
673 char *astr = mp_posit_number_tostring(mp, a_orig);
674 char *bstr = mp_posit_number_tostring(mp, b_orig);
675 snprintf(msg, 256, "Pythagorean subtraction %s+-+%s has been replaced by 0", astr, bstr);
676 mp_memory_free(astr);
677 mp_memory_free(bstr);
678 mp_error(
679 mp,
680 msg,
681 "Since I don't take square roots of negative numbers, Im zeroing this one.\n"
682 "Proceed, with fingers crossed."
683 );
684 }
685 a_orig->data.pval = mp_posit_data.zero;
686 }
687 ret->data.pval = a_orig->data.pval;
688}
689
690static void mp_posit_power_of(MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig)
691{
692 errno = 0;
693 ret->data.pval = posit_pow(a_orig->data.pval, b_orig->data.pval);
694 if (errno) {
695 mp->arith_error = 1;
696 ret->data.pval = mp_posit_data.EL_GORDO;
697 }
698}
699
700static void mp_posit_m_log(MP mp, mp_number *ret, mp_number *x_orig)
701{
702
703 if (posit_gt(x_orig->data.pval,mp_posit_data.zero)) {
704 ret->data.pval = posit_mul(posit_log(x_orig->data.pval),mp_posit_data.d256);
705 } else {
706 char msg[256];
707 char *xstr = mp_posit_number_tostring(mp, x_orig);
708 snprintf(msg, 256, "Logarithm of %s has been replaced by 0", xstr);
709 mp_memory_free(xstr);
710 mp_error(
711 mp,
712 msg,
713 "Since I don't take logs of non-positive numbers, I'm zeroing this one.\n"
714 "Proceed, with fingers crossed."
715 );
716 ret->data.pval = mp_posit_data.zero;
717 }
718}
719
720static void mp_posit_m_exp(MP mp, mp_number *ret, mp_number *x_orig)
721{
722 errno = 0;
723 ret->data.pval = posit_exp(posit_div(x_orig->data.pval,mp_posit_data.d256));
724 if (errno) {
725 if (posit_gt(x_orig->data.pval,mp_posit_data.zero)) {
726 mp->arith_error = 1;
727 ret->data.pval = mp_posit_data.EL_GORDO;
728 } else {
729 ret->data.pval = mp_posit_data.zero;
730 }
731 }
732}
733
734static void mp_posit_n_arg(MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_orig)
735{
736 if (posit_eq_zero(x_orig->data.pval) && posit_eq_zero(y_orig->data.pval)) {
737 if (posit_lt(internal_value(mp_default_zero_angle_internal).data.pval, mp_posit_data.zero)) {
738 mp_error(
739 mp,
740 "angle(0,0) is taken as zero",
741 "The 'angle' between two identical points is undefined. I'm zeroing this one.\n"
742 "Proceed, with fingers crossed."
743 );
744 ret->data.pval = mp_posit_data.zero;
745 } else {
746 ret->data.pval = internal_value(mp_default_zero_angle_internal).data.pval;
747 }
748 } else {
749 ret->type = mp_angle_type;
750
751 ret->data.pval = posit_mul(
752 posit_atan2(
753 y_orig->data.pval,
754 x_orig->data.pval
755 ),
756 mp_posit_data.d180_divided_by_pi_mul_angle
757 );
758
759 }
760}
761
762static void mp_posit_sin_cos(MP mp, mp_number *z_orig, mp_number *n_cos, mp_number *n_sin)
763{
764 posit_t rad = posit_div(z_orig->data.pval, mp_posit_data.angle_multiplier);
765 (void) mp;
766 if (posit_eq(rad, mp_posit_data.dp90) || posit_eq(rad, mp_posit_data.dm270)) {
767 n_cos->data.pval = mp_posit_data.zero;
768 n_sin->data.pval = mp_posit_data.fraction_multiplier;
769 } else if (posit_eq(rad, mp_posit_data.dm90) || posit_eq(rad, mp_posit_data.dp270)) {
770 n_cos->data.pval = mp_posit_data.zero;
771 n_sin->data.pval = mp_posit_data.negative_fraction_multiplier;
772 } else if (posit_eq(rad, mp_posit_data.dp180) || posit_eq(rad, mp_posit_data.dm180)) {
773 n_cos->data.pval = mp_posit_data.negative_fraction_multiplier;
774 n_sin->data.pval = mp_posit_data.zero;
775 } else {
776 rad = posit_mul(rad,mp_posit_data.pi_divided_by_180);
777 n_cos->data.pval = posit_mul(posit_cos(rad),mp_posit_data.fraction_multiplier);
778 n_sin->data.pval = posit_mul(posit_sin(rad),mp_posit_data.fraction_multiplier);
779 }
780}
781
782
783
784# define KK 100
785# define LL 37
786# define MM (1L<<30)
787# define mod_diff(x,y) (((x)-(y))&(MM-1))
788# define TT 70
789# define is_odd(x) ((x)&1)
790# define QUALITY 1009
791
792
793
794typedef struct mp_posit_random_info {
795 long x[KK];
796 long buf[QUALITY];
797 long dummy;
798 long started;
799 long *ptr;
800} mp_posit_random_info;
801
802static mp_posit_random_info mp_posit_random_data = {
803 .dummy = -1,
804 .started = -1,
805 .ptr = &mp_posit_random_data.dummy
806};
807
808
809
810
811static void mp_posit_aux_ran_array(long aa[], int n)
812{
813 int i, j;
814 for (j = 0; j < KK; j++) {
815 aa[j] = mp_posit_random_data.x[j];
816 }
817 for (; j < n; j++) {
818 aa[j] = mod_diff(aa[j - KK], aa[j - LL]);
819 }
820 for (i = 0; i < LL; i++, j++) {
821 mp_posit_random_data.x[i] = mod_diff(aa[j - KK], aa[j - LL]);
822 }
823 for (; i < KK; i++, j++) {
824 mp_posit_random_data.x[i] = mod_diff(aa[j - KK], mp_posit_random_data.x[i - LL]);
825 }
826}
827
828
829
830static void mp_posit_aux_ran_start(long seed)
831{
832 int t, j;
833 long x[KK + KK - 1];
834 long ss = (seed+2) & (MM - 2);
835 for (j = 0; j < KK; j++) {
836
837 x[j] = ss;
838
839 ss <<= 1;
840 if (ss >= MM) {
841 ss -= MM - 2;
842 }
843 }
844
845 x[1]++;
846 for (ss = seed & (MM - 1), t = TT - 1; t;) {
847 for (j = KK - 1; j > 0; j--) {
848
849 x[j + j] = x[j];
850 x[j + j - 1] = 0;
851 }
852 for (j = KK + KK - 2; j >= KK; j--) {
853 x[j - (KK -LL)] = mod_diff(x[j - (KK - LL)], x[j]);
854 x[j - KK] = mod_diff(x[j - KK], x[j]);
855 }
856 if (is_odd(ss)) {
857
858 for (j = KK; j>0; j--) {
859 x[j] = x[j-1];
860 }
861 x[0] = x[KK];
862
863 x[LL] = mod_diff(x[LL], x[KK]);
864 }
865 if (ss) {
866 ss >>= 1;
867 } else {
868 t--;
869 }
870 }
871 for (j = 0; j < LL; j++) {
872 mp_posit_random_data.x[j + KK - LL] = x[j];
873 }
874 for (;j < KK; j++) {
875 mp_posit_random_data.x[j - LL] = x[j];
876 }
877 for (j = 0; j < 10; j++) {
878
879 mp_posit_aux_ran_array(x, KK + KK - 1);
880 }
881 mp_posit_random_data.ptr = &mp_posit_random_data.started;
882}
883
884static long mp_posit_aux_ran_arr_cycle(void)
885{
886 if (mp_posit_random_data.ptr == &mp_posit_random_data.dummy) {
887
888 mp_posit_aux_ran_start(314159L);
889 }
890 mp_posit_aux_ran_array(mp_posit_random_data.buf, QUALITY);
891 mp_posit_random_data.buf[KK] = -1;
892 mp_posit_random_data.ptr = mp_posit_random_data.buf + 1;
893 return mp_posit_random_data.buf[0];
894}
895
896static void mp_posit_init_randoms(MP mp, int seed)
897{
898 int k = 1;
899 int j = abs(seed);
900 int f = (int) mp_fraction_multiplier;
901 while (j >= f) {
902 j = j/2;
903 }
904 for (int i = 0; i <= 54; i++) {
905 int jj = k;
906 k = j - k;
907 j = jj;
908 if (k < 0) {
909 k += f;
910 }
911 mp->randoms[(i * 21) % 55].data.pval = integer_to_posit(j);
912 }
913 mp_new_randoms(mp);
914 mp_new_randoms(mp);
915 mp_new_randoms(mp);
916
917 mp_posit_aux_ran_start((unsigned long) seed);
918}
919
920static void mp_posit_modulo(mp_number *a, mp_number *b)
921{
922 a->data.pval = posit_mul(posit_modf(posit_div(a->data.pval, b->data.pval)), b->data.pval);
923}
924
925static void mp_posit_next_unif_random(MP mp, mp_number *ret)
926{
927 unsigned long int op = (unsigned) (*mp_posit_random_data.ptr >=0 ? *mp_posit_random_data.ptr++: mp_posit_aux_ran_arr_cycle());
928 double a = op / (MM * 1.0);
929 (void) mp;
930 ret->data.pval = double_to_posit(a);
931}
932
933static void mp_posit_aux_next_random(MP mp, mp_number *ret)
934{
935 if (mp->j_random==0) {
936 mp_new_randoms(mp);
937 } else {
938 mp->j_random = mp->j_random-1;
939 }
940 mp_posit_clone(ret, &(mp->randoms[mp->j_random]));
941}
942
943static void mp_posit_m_unif_rand(MP mp, mp_number *ret, mp_number *x_orig)
944{
945 mp_number x, abs_x, u, y;
946 mp_posit_allocate_number(mp, &y, mp_fraction_type);
947 mp_posit_allocate_clone(mp, &x, mp_scaled_type, x_orig);
948 mp_posit_allocate_abs(mp, &abs_x, mp_scaled_type, &x);
949 mp_posit_allocate_number(mp, &u, mp_scaled_type);
950 mp_posit_next_unif_random(mp, &u);
951 y.data.pval = posit_mul(abs_x.data.pval, u.data.pval);
952 mp_posit_free_number(mp, &u);
953 if (mp_posit_equal(&y, &abs_x)) {
954 mp_posit_clone(ret, &((math_data *)mp->math)->md_zero_t);
955 } else if (mp_posit_greater(&x, &((math_data *)mp->math)->md_zero_t)) {
956 mp_posit_clone(ret, &y);
957 } else {
958 mp_posit_negated_clone(ret, &y);
959 }
960 mp_posit_free_number(mp, &abs_x);
961 mp_posit_free_number(mp, &x);
962 mp_posit_free_number(mp, &y);
963}
964
965static void mp_posit_m_norm_rand(MP mp, mp_number *ret)
966{
967 mp_number abs_x, u, r, la, xa;
968 mp_posit_allocate_number(mp, &la, mp_scaled_type);
969 mp_posit_allocate_number(mp, &xa, mp_scaled_type);
970 mp_posit_allocate_number(mp, &abs_x, mp_scaled_type);
971 mp_posit_allocate_number(mp, &u, mp_scaled_type);
972 mp_posit_allocate_number(mp, &r, mp_scaled_type);
973 do {
974 do {
975 mp_number v;
976 mp_posit_allocate_number(mp, &v, mp_scaled_type);
977 mp_posit_aux_next_random(mp, &v);
978 mp_posit_subtract(&v, &((math_data *)mp->math)->md_fraction_half_t);
979 mp_posit_take_fraction(mp, &xa, &((math_data *)mp->math)->md_sqrt_8_e_k, &v);
980 mp_posit_free_number(mp, &v);
981 mp_posit_aux_next_random(mp, &u);
982 mp_posit_clone(&abs_x, &xa);
983 mp_posit_abs(&abs_x);
984 } while (! mp_posit_less(&abs_x, &u));
985 mp_posit_make_fraction(mp, &r, &xa, &u);
986 mp_posit_clone(&xa, &r);
987 mp_posit_m_log(mp, &la, &u);
988 mp_posit_set_from_subtraction(&la, &((math_data *)mp->math)->md_twelve_ln_2_k, &la);
989 } while (mp_posit_ab_vs_cd(&((math_data *)mp->math)->md_one_k, &la, &xa, &xa) < 0);
990 mp_posit_clone(ret, &xa);
991 mp_posit_free_number(mp, &r);
992 mp_posit_free_number(mp, &abs_x);
993 mp_posit_free_number(mp, &la);
994 mp_posit_free_number(mp, &xa);
995 mp_posit_free_number(mp, &u);
996}
997
998static void mp_posit_set_precision(MP mp)
999{
1000 (void) mp;
1001}
1002
1003static void mp_posit_free_math(MP mp)
1004{
1005
1006 mp_posit_free_number(mp, &(mp->math->md_three_sixty_deg_t));
1007 mp_posit_free_number(mp, &(mp->math->md_one_eighty_deg_t));
1008 mp_posit_free_number(mp, &(mp->math->md_negative_one_eighty_deg_t));
1009 mp_posit_free_number(mp, &(mp->math->md_fraction_one_t));
1010 mp_posit_free_number(mp, &(mp->math->md_zero_t));
1011 mp_posit_free_number(mp, &(mp->math->md_half_unit_t));
1012 mp_posit_free_number(mp, &(mp->math->md_three_quarter_unit_t));
1013 mp_posit_free_number(mp, &(mp->math->md_unity_t));
1014 mp_posit_free_number(mp, &(mp->math->md_two_t));
1015 mp_posit_free_number(mp, &(mp->math->md_three_t));
1016 mp_posit_free_number(mp, &(mp->math->md_one_third_inf_t));
1017 mp_posit_free_number(mp, &(mp->math->md_inf_t));
1018 mp_posit_free_number(mp, &(mp->math->md_negative_inf_t));
1019 mp_posit_free_number(mp, &(mp->math->md_warning_limit_t));
1020 mp_posit_free_number(mp, &(mp->math->md_one_k));
1021 mp_posit_free_number(mp, &(mp->math->md_sqrt_8_e_k));
1022 mp_posit_free_number(mp, &(mp->math->md_twelve_ln_2_k));
1023 mp_posit_free_number(mp, &(mp->math->md_coef_bound_k));
1024 mp_posit_free_number(mp, &(mp->math->md_coef_bound_minus_1));
1025 mp_posit_free_number(mp, &(mp->math->md_fraction_threshold_t));
1026 mp_posit_free_number(mp, &(mp->math->md_half_fraction_threshold_t));
1027 mp_posit_free_number(mp, &(mp->math->md_scaled_threshold_t));
1028 mp_posit_free_number(mp, &(mp->math->md_half_scaled_threshold_t));
1029 mp_posit_free_number(mp, &(mp->math->md_near_zero_angle_t));
1030 mp_posit_free_number(mp, &(mp->math->md_p_over_v_threshold_t));
1031 mp_posit_free_number(mp, &(mp->math->md_equation_threshold_t));
1032 mp_memory_free(mp->math);
1033}
1034
1035math_data *mp_initialize_posit_math(MP mp)
1036{
1037 math_data *math = (math_data *) mp_memory_allocate(sizeof(math_data));
1038
1039 if (! mp_posit_data.initialized) {
1040 mp_posit_data.initialized = 1;
1041 mp_posit_data.unity = integer_to_posit(1);
1042 mp_posit_data.zero = integer_to_posit(0);
1043 mp_posit_data.one = integer_to_posit(1);
1044 mp_posit_data.two = integer_to_posit(2);
1045 mp_posit_data.three = integer_to_posit(3);
1046 mp_posit_data.four = integer_to_posit(4);
1047 mp_posit_data.five = integer_to_posit(5);
1048 mp_posit_data.seven = integer_to_posit(7);
1049 mp_posit_data.eight = integer_to_posit(8);
1050 mp_posit_data.sixteen = integer_to_posit(16);
1051 mp_posit_data.dp90 = integer_to_posit(90);
1052 mp_posit_data.dp180 = integer_to_posit(180);
1053 mp_posit_data.dp270 = integer_to_posit(270);
1054 mp_posit_data.dp360 = integer_to_posit(360);
1055 mp_posit_data.dm90 = integer_to_posit(-90);
1056 mp_posit_data.dm180 = integer_to_posit(-180);
1057 mp_posit_data.dm270 = integer_to_posit(-270);
1058 mp_posit_data.dm360 = integer_to_posit(-360);
1059 mp_posit_data.d16 = integer_to_posit(16);
1060 mp_posit_data.d64 = integer_to_posit(64);
1061 mp_posit_data.d256 = integer_to_posit(256);
1062 mp_posit_data.d4096 = integer_to_posit(4096);
1063 mp_posit_data.d65536 = integer_to_posit(65536);
1064 mp_posit_data.minusone = posit_neg(mp_posit_data.one);
1065 mp_posit_data.half_unit = posit_div(mp_posit_data.unity, mp_posit_data.two);
1066 mp_posit_data.three_quarter_unit = posit_mul(mp_posit_data.three, posit_div(mp_posit_data.unity,mp_posit_data.four));
1067 mp_posit_data.fraction_multiplier = integer_to_posit(mp_fraction_multiplier);
1068 mp_posit_data.negative_fraction_multiplier = posit_neg(mp_posit_data.fraction_multiplier);
1069 mp_posit_data.angle_multiplier = integer_to_posit(mp_angle_multiplier);
1070 mp_posit_data.fraction_one = mp_posit_data.fraction_multiplier;
1071 mp_posit_data.fraction_two = posit_mul(mp_posit_data.fraction_multiplier, mp_posit_data.two);
1072 mp_posit_data.fraction_three = posit_mul(mp_posit_data.fraction_multiplier, mp_posit_data.three);
1073 mp_posit_data.fraction_four = posit_mul(mp_posit_data.fraction_multiplier, mp_posit_data.four);
1074 mp_posit_data.fraction_half = posit_div(mp_posit_data.fraction_multiplier, mp_posit_data.two);
1075 mp_posit_data.fraction_one_and_half = posit_add(mp_posit_data.fraction_multiplier, mp_posit_data.fraction_half);
1076 mp_posit_data.one_eighty_degrees = posit_mul(mp_posit_data.angle_multiplier, mp_posit_data.dp180);
1077 mp_posit_data.negative_one_eighty_degrees = posit_mul(mp_posit_data.angle_multiplier, mp_posit_data.dm180);
1078 mp_posit_data.three_sixty_degrees = posit_mul(mp_posit_data.angle_multiplier, mp_posit_data.dp360);
1079 mp_posit_data.no_crossing = posit_add(mp_posit_data.fraction_multiplier, mp_posit_data.one);
1080 mp_posit_data.one_crossing = mp_posit_data.fraction_multiplier;
1081 mp_posit_data.zero_crossing = mp_posit_data.zero;
1082 mp_posit_data.error_correction = double_to_posit(1E-12);
1083 mp_posit_data.warning_limit = posit_pow(mp_posit_data.two, integer_to_posit(52));
1084 mp_posit_data.pi = double_to_posit(3.1415926535897932384626433832795028841971);
1085 mp_posit_data.pi_divided_by_180 = posit_div(mp_posit_data.pi, mp_posit_data.dp180);
1086 mp_posit_data.epsilon = posit_pow(mp_posit_data.two, integer_to_posit(-52.0));
1087 mp_posit_data.EL_GORDO = posit_sub(posit_div(double_to_posit(DBL_MAX),mp_posit_data.two), mp_posit_data.one);
1088 mp_posit_data.negative_EL_GORDO = posit_neg(mp_posit_data.EL_GORDO);
1089 mp_posit_data.one_third_EL_GORDO = posit_div(mp_posit_data.EL_GORDO, mp_posit_data.three);
1090 mp_posit_data.coef = posit_div(mp_posit_data.seven, mp_posit_data.three);
1091 mp_posit_data.coef_bound = posit_mul(mp_posit_data.coef, mp_posit_data.fraction_multiplier);
1092 mp_posit_data.scaled_threshold = double_to_posit(0.000122);
1093 mp_posit_data.near_zero_angle = posit_mul(double_to_posit(0.0256), mp_posit_data.angle_multiplier);
1094 mp_posit_data.p_over_v_threshold = integer_to_posit(0x80000);
1095 mp_posit_data.equation_threshold = double_to_posit(0.001);
1096 mp_posit_data.sqrt_two_mul_fraction_one =
1097 posit_mul(
1098 posit_sqrt(mp_posit_data.two),
1099 mp_posit_data.fraction_one
1100 );
1101 mp_posit_data.sqrt_five_minus_one_mul_fraction_one_and_half =
1102 posit_mul(
1103 posit_mul(
1104 mp_posit_data.three,
1105 mp_posit_data.fraction_half
1106 ),
1107 posit_sub(
1108 posit_sqrt(mp_posit_data.five),
1109 mp_posit_data.one
1110 )
1111 );
1112 mp_posit_data.three_minus_sqrt_five_mul_fraction_one_and_half =
1113 posit_mul(
1114 posit_mul(
1115 mp_posit_data.three,
1116 mp_posit_data.fraction_half
1117 ),
1118 posit_sub(
1119 mp_posit_data.three,
1120 posit_sqrt(mp_posit_data.five)
1121 )
1122 );
1123 mp_posit_data.d180_divided_by_pi_mul_angle =
1124 posit_mul(
1125 posit_div(
1126 mp_posit_data.dp180,
1127 mp_posit_data.pi
1128 ),
1129 mp_posit_data.angle_multiplier
1130 );
1131 }
1132
1133 math->md_allocate = mp_posit_allocate_number;
1134 math->md_free = mp_posit_free_number;
1135 math->md_allocate_clone = mp_posit_allocate_clone;
1136 math->md_allocate_abs = mp_posit_allocate_abs;
1137 math->md_allocate_double = mp_posit_allocate_double;
1138
1139 mp_posit_allocate_number(mp, &math->md_precision_default, mp_scaled_type);
1140 mp_posit_allocate_number(mp, &math->md_precision_max, mp_scaled_type);
1141 mp_posit_allocate_number(mp, &math->md_precision_min, mp_scaled_type);
1142
1143 mp_posit_allocate_number(mp, &math->md_epsilon_t, mp_scaled_type);
1144 mp_posit_allocate_number(mp, &math->md_inf_t, mp_scaled_type);
1145 mp_posit_allocate_number(mp, &math->md_negative_inf_t, mp_scaled_type);
1146 mp_posit_allocate_number(mp, &math->md_warning_limit_t, mp_scaled_type);
1147 mp_posit_allocate_number(mp, &math->md_one_third_inf_t, mp_scaled_type);
1148 mp_posit_allocate_number(mp, &math->md_unity_t, mp_scaled_type);
1149 mp_posit_allocate_number(mp, &math->md_two_t, mp_scaled_type);
1150 mp_posit_allocate_number(mp, &math->md_three_t, mp_scaled_type);
1151 mp_posit_allocate_number(mp, &math->md_half_unit_t, mp_scaled_type);
1152 mp_posit_allocate_number(mp, &math->md_three_quarter_unit_t, mp_scaled_type);
1153 mp_posit_allocate_number(mp, &math->md_zero_t, mp_scaled_type);
1154
1155 mp_posit_allocate_number(mp, &math->md_arc_tol_k, mp_fraction_type);
1156 mp_posit_allocate_number(mp, &math->md_fraction_one_t, mp_fraction_type);
1157 mp_posit_allocate_number(mp, &math->md_fraction_half_t, mp_fraction_type);
1158 mp_posit_allocate_number(mp, &math->md_fraction_three_t, mp_fraction_type);
1159 mp_posit_allocate_number(mp, &math->md_fraction_four_t, mp_fraction_type);
1160
1161 mp_posit_allocate_number(mp, &math->md_three_sixty_deg_t, mp_angle_type);
1162 mp_posit_allocate_number(mp, &math->md_one_eighty_deg_t, mp_angle_type);
1163 mp_posit_allocate_number(mp, &math->md_negative_one_eighty_deg_t, mp_angle_type);
1164
1165 mp_posit_allocate_number(mp, &math->md_one_k, mp_scaled_type);
1166 mp_posit_allocate_number(mp, &math->md_sqrt_8_e_k, mp_scaled_type);
1167 mp_posit_allocate_number(mp, &math->md_twelve_ln_2_k, mp_fraction_type);
1168 mp_posit_allocate_number(mp, &math->md_coef_bound_k, mp_fraction_type);
1169 mp_posit_allocate_number(mp, &math->md_coef_bound_minus_1, mp_fraction_type);
1170 mp_posit_allocate_number(mp, &math->md_twelvebits_3, mp_scaled_type);
1171 mp_posit_allocate_number(mp, &math->md_twentysixbits_sqrt2_t, mp_fraction_type);
1172 mp_posit_allocate_number(mp, &math->md_twentyeightbits_d_t, mp_fraction_type);
1173 mp_posit_allocate_number(mp, &math->md_twentysevenbits_sqrt2_d_t, mp_fraction_type);
1174
1175 mp_posit_allocate_number(mp, &math->md_fraction_threshold_t, mp_fraction_type);
1176 mp_posit_allocate_number(mp, &math->md_half_fraction_threshold_t, mp_fraction_type);
1177 mp_posit_allocate_number(mp, &math->md_scaled_threshold_t, mp_scaled_type);
1178 mp_posit_allocate_number(mp, &math->md_half_scaled_threshold_t, mp_scaled_type);
1179 mp_posit_allocate_number(mp, &math->md_near_zero_angle_t, mp_angle_type);
1180 mp_posit_allocate_number(mp, &math->md_p_over_v_threshold_t, mp_fraction_type);
1181 mp_posit_allocate_number(mp, &math->md_equation_threshold_t, mp_scaled_type);
1182
1183 math->md_precision_default.data.pval = posit_mul(mp_posit_data.d16, mp_posit_data.unity);
1184 math->md_precision_max.data.pval = posit_mul(mp_posit_data.d16, mp_posit_data.unity);
1185 math->md_precision_min.data.pval = posit_mul(mp_posit_data.d16, mp_posit_data.unity);
1186 math->md_epsilon_t.data.pval = mp_posit_data.epsilon;
1187 math->md_inf_t.data.pval = mp_posit_data.EL_GORDO;
1188 math->md_negative_inf_t.data.pval = mp_posit_data.negative_EL_GORDO;
1189 math->md_one_third_inf_t.data.pval = mp_posit_data.one_third_EL_GORDO;
1190 math->md_warning_limit_t.data.pval = mp_posit_data.warning_limit;
1191 math->md_unity_t.data.pval = mp_posit_data.unity;
1192 math->md_two_t.data.pval = mp_posit_data.two;
1193 math->md_three_t.data.pval = mp_posit_data.three;
1194 math->md_half_unit_t.data.pval = mp_posit_data.half_unit;
1195 math->md_three_quarter_unit_t.data.pval = mp_posit_data.three_quarter_unit;
1196 math->md_arc_tol_k.data.pval = posit_div(mp_posit_data.unity, mp_posit_data.d4096);
1197 math->md_fraction_one_t.data.pval = mp_posit_data.fraction_one;
1198 math->md_fraction_half_t.data.pval = mp_posit_data.fraction_half;
1199 math->md_fraction_three_t.data.pval = mp_posit_data.fraction_three;
1200 math->md_fraction_four_t.data.pval = mp_posit_data.fraction_four;
1201 math->md_three_sixty_deg_t.data.pval = mp_posit_data.three_sixty_degrees;
1202 math->md_one_eighty_deg_t.data.pval = mp_posit_data.one_eighty_degrees;
1203 math->md_negative_one_eighty_deg_t.data.pval = mp_posit_data.negative_one_eighty_degrees;
1204 math->md_one_k.data.pval = posit_div(mp_posit_data.one, mp_posit_data.d64);
1205 math->md_sqrt_8_e_k.data.pval = double_to_posit(1.71552776992141359295);
1206 math->md_twelve_ln_2_k.data.pval = posit_mul(double_to_posit(8.31776616671934371292), mp_posit_data.d256);
1207 math->md_twelvebits_3.data.pval = posit_div(integer_to_posit(1365), mp_posit_data.unity);
1208 math->md_twentysixbits_sqrt2_t.data.pval = posit_div(integer_to_posit(94906266), mp_posit_data.d65536);
1209 math->md_twentyeightbits_d_t.data.pval = posit_div(integer_to_posit(35596755), mp_posit_data.d65536);
1210 math->md_twentysevenbits_sqrt2_d_t.data.pval = posit_div(integer_to_posit(25170707), mp_posit_data.d65536);
1211 math->md_coef_bound_k.data.pval = mp_posit_data.coef_bound;
1212 math->md_coef_bound_minus_1.data.pval = posit_sub(mp_posit_data.coef_bound, posit_div(mp_posit_data.one, mp_posit_data.d65536));
1213 math->md_fraction_threshold_t.data.pval = double_to_posit(0.04096);
1214 math->md_half_fraction_threshold_t.data.pval = posit_div(mp_posit_data.fraction_threshold, mp_posit_data.two);
1215 math->md_scaled_threshold_t.data.pval = mp_posit_data.scaled_threshold;
1216 math->md_half_scaled_threshold_t.data.pval = posit_div(mp_posit_data.scaled_threshold,mp_posit_data.two);
1217 math->md_near_zero_angle_t.data.pval = mp_posit_data.near_zero_angle;
1218 math->md_p_over_v_threshold_t.data.pval = mp_posit_data.p_over_v_threshold;
1219 math->md_equation_threshold_t.data.pval = mp_posit_data.equation_threshold;
1220
1221 math->md_from_int = mp_posit_set_from_int;
1222 math->md_from_boolean = mp_posit_set_from_boolean;
1223 math->md_from_scaled = mp_posit_set_from_scaled;
1224 math->md_from_double = mp_posit_set_from_double;
1225 math->md_from_addition = mp_posit_set_from_addition;
1226 math->md_half_from_addition = mp_posit_set_half_from_addition;
1227 math->md_from_subtraction = mp_posit_set_from_subtraction;
1228 math->md_half_from_subtraction = mp_posit_set_half_from_subtraction;
1229 math->md_from_of_the_way = mp_posit_set_from_of_the_way;
1230 math->md_from_div = mp_posit_set_from_div;
1231 math->md_from_mul = mp_posit_set_from_mul;
1232 math->md_from_int_div = mp_posit_set_from_int_div;
1233 math->md_from_int_mul = mp_posit_set_from_int_mul;
1234 math->md_negate = mp_posit_negate;
1235 math->md_add = mp_posit_add;
1236 math->md_subtract = mp_posit_subtract;
1237 math->md_half = mp_posit_half;
1238 math->md_do_double = mp_posit_double;
1239 math->md_abs = mp_posit_abs;
1240 math->md_clone = mp_posit_clone;
1241 math->md_negated_clone = mp_posit_negated_clone;
1242 math->md_abs_clone = mp_posit_abs_clone;
1243 math->md_swap = mp_posit_swap;
1244 math->md_add_scaled = mp_posit_add_scaled;
1245 math->md_multiply_int = mp_posit_multiply_int;
1246 math->md_divide_int = mp_posit_divide_int;
1247 math->md_to_int = mp_posit_to_int;
1248 math->md_to_boolean = mp_posit_to_boolean;
1249 math->md_to_scaled = mp_posit_to_scaled;
1250 math->md_to_double = mp_posit_to_double;
1251 math->md_odd = mp_posit_odd;
1252 math->md_equal = mp_posit_equal;
1253 math->md_less = mp_posit_less;
1254 math->md_greater = mp_posit_greater;
1255 math->md_non_equal_abs = mp_posit_non_equal_abs;
1256 math->md_round_unscaled = mp_posit_round_unscaled;
1257 math->md_floor_scaled = mp_posit_floor;
1258 math->md_fraction_to_round_scaled = mp_posit_fraction_to_round_scaled;
1259 math->md_make_scaled = mp_posit_make_scaled;
1260 math->md_make_fraction = mp_posit_make_fraction;
1261 math->md_take_fraction = mp_posit_take_fraction;
1262 math->md_take_scaled = mp_posit_take_scaled;
1263 math->md_velocity = mp_posit_velocity;
1264 math->md_n_arg = mp_posit_n_arg;
1265 math->md_m_log = mp_posit_m_log;
1266 math->md_m_exp = mp_posit_m_exp;
1267 math->md_m_unif_rand = mp_posit_m_unif_rand;
1268 math->md_m_norm_rand = mp_posit_m_norm_rand;
1269 math->md_pyth_add = mp_posit_pyth_add;
1270 math->md_pyth_sub = mp_posit_pyth_sub;
1271 math->md_power_of = mp_posit_power_of;
1272 math->md_fraction_to_scaled = mp_posit_fraction_to_scaled;
1273 math->md_scaled_to_fraction = mp_posit_scaled_to_fraction;
1274 math->md_scaled_to_angle = mp_posit_scaled_to_angle;
1275 math->md_angle_to_scaled = mp_posit_angle_to_scaled;
1276 math->md_init_randoms = mp_posit_init_randoms;
1277 math->md_sin_cos = mp_posit_sin_cos;
1278 math->md_slow_add = mp_posit_slow_add;
1279 math->md_sqrt = mp_posit_square_rt;
1280 math->md_print = mp_posit_print_number;
1281 math->md_tostring = mp_posit_number_tostring;
1282 math->md_modulo = mp_posit_modulo;
1283 math->md_ab_vs_cd = mp_posit_ab_vs_cd;
1284 math->md_crossing_point = mp_posit_crossing_point;
1285 math->md_scan_numeric = mp_posit_scan_numeric_token;
1286 math->md_scan_fractional = mp_posit_scan_fractional_token;
1287
1288 math->md_free_math = mp_posit_free_math;
1289 math->md_set_precision = mp_posit_set_precision;
1290 return math;
1291}
1292 |