1 /* Copyright (C) 2001-2019 Artifex Software, Inc.
2 All Rights Reserved.
3
4 This software is provided AS-IS with no warranty, either express or
5 implied.
6
7 This software is distributed under license and may not be copied,
8 modified or distributed except as expressly authorized under the terms
9 of the license contained in the file LICENSE in this distribution.
10
11 Refer to licensing information at http://www.artifex.com or contact
12 Artifex Software, Inc., 1305 Grant Avenue - Suite 200, Novato,
13 CA 94945, U.S.A., +1(415)492-9861, for further information.
14 */
15
16
17 /* Double-precision floating point arithmetic operators */
18 #include "math_.h"
19 #include "memory_.h"
20 #include "string_.h"
21 #include "ctype_.h"
22 #include "ghost.h"
23 #include "gxfarith.h"
24 #include "oper.h"
25 #include "store.h"
26
27 /*
28 * Thanks to Jean-Pierre Demailly of the Institut Fourier of the
29 * Universit\'e de Grenoble I <demailly@fourier.grenet.fr> for proposing
30 * this package and for arranging the funding for its creation.
31 *
32 * These operators work with doubles represented as 8-byte strings. When
33 * applicable, they write their result into a string supplied as an argument.
34 * They also accept ints and reals as arguments.
35 */
36
37 /* Forward references */
38 static int double_params_result(os_ptr, int, double *);
39 static int double_params(os_ptr, int, double *);
40 static int double_result(i_ctx_t *, int, double);
41 static int double_unary(i_ctx_t *, double (*)(double));
42
43 #define dbegin_unary()\
44 os_ptr op = osp;\
45 double num;\
46 int code = double_params_result(op, 1, &num);\
47 \
48 if ( code < 0 )\
49 return code
50
51 #define dbegin_binary()\
52 os_ptr op = osp;\
53 double num[2];\
54 int code = double_params_result(op, 2, num);\
55 \
56 if ( code < 0 )\
57 return code
58
59 /* ------ Arithmetic ------ */
60
61 /* <dnum1> <dnum2> <dresult> .dadd <dresult> */
62 static int
zdadd(i_ctx_t * i_ctx_p)63 zdadd(i_ctx_t *i_ctx_p)
64 {
65 dbegin_binary();
66 return double_result(i_ctx_p, 2, num[0] + num[1]);
67 }
68
69 /* <dnum1> <dnum2> <dresult> .ddiv <dresult> */
70 static int
zddiv(i_ctx_t * i_ctx_p)71 zddiv(i_ctx_t *i_ctx_p)
72 {
73 dbegin_binary();
74 if (num[1] == 0.0)
75 return_error(gs_error_undefinedresult);
76 return double_result(i_ctx_p, 2, num[0] / num[1]);
77 }
78
79 /* <dnum1> <dnum2> <dresult> .dmul <dresult> */
80 static int
zdmul(i_ctx_t * i_ctx_p)81 zdmul(i_ctx_t *i_ctx_p)
82 {
83 dbegin_binary();
84 return double_result(i_ctx_p, 2, num[0] * num[1]);
85 }
86
87 /* <dnum1> <dnum2> <dresult> .dsub <dresult> */
88 static int
zdsub(i_ctx_t * i_ctx_p)89 zdsub(i_ctx_t *i_ctx_p)
90 {
91 dbegin_binary();
92 return double_result(i_ctx_p, 2, num[0] - num[1]);
93 }
94
95 /* ------ Simple functions ------ */
96
97 /* <dnum> <dresult> .dabs <dresult> */
98 static int
zdabs(i_ctx_t * i_ctx_p)99 zdabs(i_ctx_t *i_ctx_p)
100 {
101 return double_unary(i_ctx_p, fabs);
102 }
103
104 /* <dnum> <dresult> .dceiling <dresult> */
105 static int
zdceiling(i_ctx_t * i_ctx_p)106 zdceiling(i_ctx_t *i_ctx_p)
107 {
108 return double_unary(i_ctx_p, ceil);
109 }
110
111 /* <dnum> <dresult> .dfloor <dresult> */
112 static int
zdfloor(i_ctx_t * i_ctx_p)113 zdfloor(i_ctx_t *i_ctx_p)
114 {
115 return double_unary(i_ctx_p, floor);
116 }
117
118 /* <dnum> <dresult> .dneg <dresult> */
119 static int
zdneg(i_ctx_t * i_ctx_p)120 zdneg(i_ctx_t *i_ctx_p)
121 {
122 dbegin_unary();
123 return double_result(i_ctx_p, 1, -num);
124 }
125
126 /* <dnum> <dresult> .dround <dresult> */
127 static int
zdround(i_ctx_t * i_ctx_p)128 zdround(i_ctx_t *i_ctx_p)
129 {
130 dbegin_unary();
131 return double_result(i_ctx_p, 1, floor(num + 0.5));
132 }
133
134 /* <dnum> <dresult> .dsqrt <dresult> */
135 static int
zdsqrt(i_ctx_t * i_ctx_p)136 zdsqrt(i_ctx_t *i_ctx_p)
137 {
138 dbegin_unary();
139 if (num < 0.0)
140 return_error(gs_error_rangecheck);
141 return double_result(i_ctx_p, 1, sqrt(num));
142 }
143
144 /* <dnum> <dresult> .dtruncate <dresult> */
145 static int
zdtruncate(i_ctx_t * i_ctx_p)146 zdtruncate(i_ctx_t *i_ctx_p)
147 {
148 dbegin_unary();
149 return double_result(i_ctx_p, 1, (num < 0 ? ceil(num) : floor(num)));
150 }
151
152 /* ------ Transcendental functions ------ */
153
154 static int
darc(i_ctx_t * i_ctx_p,double (* afunc)(double))155 darc(i_ctx_t *i_ctx_p, double (*afunc)(double))
156 {
157 dbegin_unary();
158 return double_result(i_ctx_p, 1, (*afunc)(num) * radians_to_degrees);
159 }
160 /* <dnum> <dresult> .darccos <dresult> */
161 static int
zdarccos(i_ctx_t * i_ctx_p)162 zdarccos(i_ctx_t *i_ctx_p)
163 {
164 return darc(i_ctx_p, acos);
165 }
166 /* <dnum> <dresult> .darcsin <dresult> */
167 static int
zdarcsin(i_ctx_t * i_ctx_p)168 zdarcsin(i_ctx_t *i_ctx_p)
169 {
170 return darc(i_ctx_p, asin);
171 }
172
173 /* <dnum> <ddenom> <dresult> .datan <dresult> */
174 static int
zdatan(i_ctx_t * i_ctx_p)175 zdatan(i_ctx_t *i_ctx_p)
176 {
177 double result;
178
179 dbegin_binary();
180 if (num[0] == 0) { /* on X-axis, special case */
181 if (num[1] == 0)
182 return_error(gs_error_undefinedresult);
183 result = (num[1] < 0 ? 180 : 0);
184 } else {
185 result = atan2(num[0], num[1]) * radians_to_degrees;
186 if (result < 0)
187 result += 360;
188 }
189 return double_result(i_ctx_p, 2, result);
190 }
191
192 /* <dnum> <dresult> .dcos <dresult> */
193 static int
zdcos(i_ctx_t * i_ctx_p)194 zdcos(i_ctx_t *i_ctx_p)
195 {
196 return double_unary(i_ctx_p, gs_cos_degrees);
197 }
198
199 /* <dbase> <dexponent> <dresult> .dexp <dresult> */
200 static int
zdexp(i_ctx_t * i_ctx_p)201 zdexp(i_ctx_t *i_ctx_p)
202 {
203 double ipart;
204
205 dbegin_binary();
206 if (num[0] == 0.0 && num[1] == 0.0)
207 return_error(gs_error_undefinedresult);
208 if (num[0] < 0.0 && modf(num[1], &ipart) != 0.0)
209 return_error(gs_error_undefinedresult);
210 return double_result(i_ctx_p, 2, pow(num[0], num[1]));
211 }
212
213 static int
dlog(i_ctx_t * i_ctx_p,double (* lfunc)(double))214 dlog(i_ctx_t *i_ctx_p, double (*lfunc)(double))
215 {
216 dbegin_unary();
217 if (num <= 0.0)
218 return_error(gs_error_rangecheck);
219 return double_result(i_ctx_p, 1, (*lfunc)(num));
220 }
221 /* <dposnum> <dresult> .dln <dresult> */
222 static int
zdln(i_ctx_t * i_ctx_p)223 zdln(i_ctx_t *i_ctx_p)
224 {
225 return dlog(i_ctx_p, log);
226 }
227 /* <dposnum> <dresult> .dlog <dresult> */
228 static int
zdlog(i_ctx_t * i_ctx_p)229 zdlog(i_ctx_t *i_ctx_p)
230 {
231 return dlog(i_ctx_p, log10);
232 }
233
234 /* <dnum> <dresult> .dsin <dresult> */
235 static int
zdsin(i_ctx_t * i_ctx_p)236 zdsin(i_ctx_t *i_ctx_p)
237 {
238 return double_unary(i_ctx_p, gs_sin_degrees);
239 }
240
241 /* ------ Comparison ------ */
242
243 static int
dcompare(i_ctx_t * i_ctx_p,int mask)244 dcompare(i_ctx_t *i_ctx_p, int mask)
245 {
246 os_ptr op = osp;
247 double num[2];
248 int code = double_params(op, 2, num);
249
250 if (code < 0)
251 return code;
252 make_bool(op - 1,
253 (mask & (num[0] < num[1] ? 1 : num[0] > num[1] ? 4 : 2))
254 != 0);
255 pop(1);
256 return 0;
257 }
258 /* <dnum1> <dnum2> .deq <bool> */
259 static int
zdeq(i_ctx_t * i_ctx_p)260 zdeq(i_ctx_t *i_ctx_p)
261 {
262 return dcompare(i_ctx_p, 2);
263 }
264 /* <dnum1> <dnum2> .dge <bool> */
265 static int
zdge(i_ctx_t * i_ctx_p)266 zdge(i_ctx_t *i_ctx_p)
267 {
268 return dcompare(i_ctx_p, 6);
269 }
270 /* <dnum1> <dnum2> .dgt <bool> */
271 static int
zdgt(i_ctx_t * i_ctx_p)272 zdgt(i_ctx_t *i_ctx_p)
273 {
274 return dcompare(i_ctx_p, 4);
275 }
276 /* <dnum1> <dnum2> .dle <bool> */
277 static int
zdle(i_ctx_t * i_ctx_p)278 zdle(i_ctx_t *i_ctx_p)
279 {
280 return dcompare(i_ctx_p, 3);
281 }
282 /* <dnum1> <dnum2> .dlt <bool> */
283 static int
zdlt(i_ctx_t * i_ctx_p)284 zdlt(i_ctx_t *i_ctx_p)
285 {
286 return dcompare(i_ctx_p, 1);
287 }
288 /* <dnum1> <dnum2> .dne <bool> */
289 static int
zdne(i_ctx_t * i_ctx_p)290 zdne(i_ctx_t *i_ctx_p)
291 {
292 return dcompare(i_ctx_p, 5);
293 }
294
295 /* ------ Conversion ------ */
296
297 /* Take the easy way out.... */
298 #define MAX_CHARS 50
299
300 /* <dnum> <dresult> .cvd <dresult> */
301 static int
zcvd(i_ctx_t * i_ctx_p)302 zcvd(i_ctx_t *i_ctx_p)
303 {
304 dbegin_unary();
305 return double_result(i_ctx_p, 1, num);
306 }
307
308 /* <string> <dresult> .cvsd <dresult> */
309 static int
zcvsd(i_ctx_t * i_ctx_p)310 zcvsd(i_ctx_t *i_ctx_p)
311 {
312 os_ptr op = osp;
313 int code = double_params_result(op, 0, NULL);
314 double num;
315 char dot, buf[MAX_CHARS + 2];
316 char *str = buf;
317 uint len;
318 char end;
319
320 if (code < 0)
321 return code;
322 check_read_type(op[-1], t_string);
323 len = r_size(op - 1);
324 if (len > MAX_CHARS)
325 return_error(gs_error_limitcheck);
326 gs_sprintf(buf, "%f", 1.5);
327 dot = buf[1]; /* locale-dependent */
328 memcpy(str, op[-1].value.bytes, len);
329 /*
330 * We check syntax in the following way: we remove whitespace,
331 * verify that the string contains only [0123456789+-.dDeE],
332 * then append a $ and then check that the next character after
333 * the scanned number is a $.
334 */
335 while (len > 0 && isspace(*str))
336 ++str, --len;
337 while (len > 0 && isspace(str[len - 1]))
338 --len;
339 str[len] = 0;
340 if (strspn(str, "0123456789+-.dDeE") != len)
341 return_error(gs_error_syntaxerror);
342 strcat(str, "$");
343 if (dot != '.') {
344 char *pdot = strchr(str, '.');
345 if (pdot)
346 *pdot = dot;
347 }
348 if (sscanf(str, "%lf%c", &num, &end) != 2 || end != '$')
349 return_error(gs_error_syntaxerror);
350 return double_result(i_ctx_p, 1, num);
351 }
352
353 /* <dnum> .dcvi <int> */
354 static int
zdcvi(i_ctx_t * i_ctx_p)355 zdcvi(i_ctx_t *i_ctx_p)
356 {
357 os_ptr op = osp;
358 #define alt_min_long (-1L << (ARCH_SIZEOF_LONG * 8 - 1))
359 #define alt_max_long (~(alt_min_long))
360 static const double min_int_real = (alt_min_long * 1.0 - 1);
361 static const double max_int_real = (alt_max_long * 1.0 + 1);
362 double num;
363 int code = double_params(op, 1, &num);
364
365 if (code < 0)
366 return code;
367
368 if (num < min_int_real || num > max_int_real)
369 return_error(gs_error_rangecheck);
370 make_int(op, (long)num); /* truncates toward 0 */
371 return 0;
372 }
373
374 /* <dnum> .dcvr <real> */
375 static int
zdcvr(i_ctx_t * i_ctx_p)376 zdcvr(i_ctx_t *i_ctx_p)
377 {
378 os_ptr op = osp;
379 #define b30 (0x40000000L * 1.0)
380 #define max_mag (0xffffff * b30 * b30 * b30 * 0x4000)
381 static const float min_real = -max_mag;
382 static const float max_real = max_mag;
383 #undef b30
384 #undef max_mag
385 double num;
386 int code = double_params(op, 1, &num);
387
388 if (code < 0)
389 return code;
390 if (num < min_real || num > max_real)
391 return_error(gs_error_rangecheck);
392 make_real(op, (float)num);
393 return 0;
394 }
395
396 /* <dnum> <string> .dcvs <substring> */
397 static int
zdcvs(i_ctx_t * i_ctx_p)398 zdcvs(i_ctx_t *i_ctx_p)
399 {
400 os_ptr op = osp;
401 double num;
402 int code = double_params(op - 1, 1, &num);
403 char dot, str[MAX_CHARS + 1];
404 int len;
405
406 if (code < 0)
407 return code;
408 check_write_type(*op, t_string);
409 gs_sprintf(str, "%f", 1.5);
410 dot = str[1]; /* locale-dependent */
411 /*
412 * To get fully accurate output results for IEEE double-
413 * precision floats (53 bits of mantissa), the ANSI
414 * %g default of 6 digits is not enough; 16 are needed.
415 * Unfortunately, using %.16g produces unfortunate artifacts such as
416 * 1.2 printing as 1.200000000000005. Therefore, we print using %g,
417 * and if the result isn't accurate enough, print again
418 * using %.16g.
419 */
420 {
421 double scanned;
422
423 gs_sprintf(str, "%g", num);
424 sscanf(str, "%lf", &scanned);
425 if (scanned != num)
426 gs_sprintf(str, "%.16g", num);
427 }
428 len = strlen(str);
429 if (len > r_size(op))
430 return_error(gs_error_rangecheck);
431 /* Juggling locales isn't thread-safe. Posix me harder. */
432 if (dot != '.') {
433 char *pdot = strchr(str, dot);
434 if (pdot)
435 *pdot = '.';
436 }
437 memcpy(op->value.bytes, str, len);
438 op[-1] = *op;
439 r_set_size(op - 1, len);
440 pop(1);
441 return 0;
442 }
443
444 /* ------ Initialization table ------ */
445
446 /* We need to split the table because of the 16-element limit. */
447 const op_def zdouble1_op_defs[] = {
448 /* Arithmetic */
449 {"3.dadd", zdadd},
450 {"3.ddiv", zddiv},
451 {"3.dmul", zdmul},
452 {"3.dsub", zdsub},
453 /* Comparison */
454 {"2.deq", zdeq},
455 {"2.dge", zdge},
456 {"2.dgt", zdgt},
457 {"2.dle", zdle},
458 {"2.dlt", zdlt},
459 {"2.dne", zdne},
460 /* Conversion */
461 {"2.cvd", zcvd},
462 {"2.cvsd", zcvsd},
463 {"1.dcvi", zdcvi},
464 {"1.dcvr", zdcvr},
465 {"2.dcvs", zdcvs},
466 op_def_end(0)
467 };
468 const op_def zdouble2_op_defs[] = {
469 /* Simple functions */
470 {"2.dabs", zdabs},
471 {"2.dceiling", zdceiling},
472 {"2.dfloor", zdfloor},
473 {"2.dneg", zdneg},
474 {"2.dround", zdround},
475 {"2.dsqrt", zdsqrt},
476 {"2.dtruncate", zdtruncate},
477 /* Transcendental functions */
478 {"2.darccos", zdarccos},
479 {"2.darcsin", zdarcsin},
480 {"3.datan", zdatan},
481 {"2.dcos", zdcos},
482 {"3.dexp", zdexp},
483 {"2.dln", zdln},
484 {"2.dlog", zdlog},
485 {"2.dsin", zdsin},
486 op_def_end(0)
487 };
488
489 /* ------ Internal procedures ------ */
490
491 /* Get some double arguments. */
492 static int
double_params(os_ptr op,int count,double * pval)493 double_params(os_ptr op, int count, double *pval)
494 {
495 pval += count;
496 while (--count >= 0) {
497 switch (r_type(op)) {
498 case t_real:
499 *--pval = op->value.realval;
500 break;
501 case t_integer:
502 *--pval = op->value.intval;
503 break;
504 case t_string:
505 if (!r_has_attr(op, a_read) ||
506 r_size(op) != sizeof(double)
507 )
508 return_error(gs_error_typecheck);
509 --pval;
510 memcpy(pval, op->value.bytes, sizeof(double));
511 break;
512 case t__invalid:
513 return_error(gs_error_stackunderflow);
514 default:
515 return_error(gs_error_typecheck);
516 }
517 op--;
518 }
519 return 0;
520 }
521
522 /* Get some double arguments, and check for a double result. */
523 static int
double_params_result(os_ptr op,int count,double * pval)524 double_params_result(os_ptr op, int count, double *pval)
525 {
526 check_write_type(*op, t_string);
527 if (r_size(op) != sizeof(double))
528 return_error(gs_error_typecheck);
529 return double_params(op - 1, count, pval);
530 }
531
532 /* Return a double result. */
533 static int
double_result(i_ctx_t * i_ctx_p,int count,double result)534 double_result(i_ctx_t *i_ctx_p, int count, double result)
535 {
536 os_ptr op = osp;
537 os_ptr op1 = op - count;
538
539 ref_assign_inline(op1, op);
540 memcpy(op1->value.bytes, &result, sizeof(double));
541 pop(count);
542 return 0;
543 }
544
545 /* Apply a unary function to a double operand. */
546 static int
double_unary(i_ctx_t * i_ctx_p,double (* func)(double))547 double_unary(i_ctx_t *i_ctx_p, double (*func)(double))
548 {
549 dbegin_unary();
550 return double_result(i_ctx_p, 1, (*func)(num));
551 }
552