1 /* Copyright (C) 1989, 2000 Aladdin Enterprises. All rights reserved.
2
3 This software is provided AS-IS with no warranty, either express or
4 implied.
5
6 This software is distributed under license and may not be copied,
7 modified or distributed except as expressly authorized under the terms
8 of the license contained in the file LICENSE in this distribution.
9
10 For more information about licensing, please refer to
11 http://www.ghostscript.com/licensing/. For information on
12 commercial licensing, go to http://www.artifex.com/licensing/ or
13 contact Artifex Software, Inc., 101 Lucas Valley Road #110,
14 San Rafael, CA 94903, U.S.A., +1(415)492-9861.
15 */
16
17 /* $Id: gsmisc.c,v 1.22 2004/12/08 21:35:13 stefan Exp $ */
18 /* Miscellaneous utilities for Ghostscript library */
19
20 /*
21 * In order to capture the original definition of sqrt, which might be
22 * either a procedure or a macro and might not have an ANSI-compliant
23 * prototype (!), we need to do the following:
24 */
25 #include "std.h"
26 #if defined(VMS) && defined(__GNUC__)
27 /* DEC VAX/VMS C comes with a math.h file, but GNU VAX/VMS C does not. */
28 # include "vmsmath.h"
29 #else
30 # include <math.h>
31 #endif
32 inline private double
orig_sqrt(double x)33 orig_sqrt(double x)
34 {
35 return sqrt(x);
36 }
37
38 /* Here is the real #include section. */
39 #include "ctype_.h"
40 #include "malloc_.h"
41 #include "math_.h"
42 #include "memory_.h"
43 #include "string_.h"
44 #include "gx.h"
45 #include "gpcheck.h" /* for gs_return_check_interrupt */
46 #include "gserror.h" /* for prototype */
47 #include "gserrors.h"
48 #include "gconfigv.h" /* for USE_ASM */
49 #include "gxfarith.h"
50 #include "gxfixed.h"
51
52 /* ------ Redirected stdout and stderr ------ */
53
54 #include <stdarg.h>
55 #define PRINTF_BUF_LENGTH 1024
56
outprintf(const gs_memory_t * mem,const char * fmt,...)57 int outprintf(const gs_memory_t *mem, const char *fmt, ...)
58 {
59 int count;
60 char buf[PRINTF_BUF_LENGTH];
61 va_list args;
62
63 va_start(args, fmt);
64
65 count = vsprintf(buf, fmt, args);
66 outwrite(mem, buf, count);
67 if (count >= PRINTF_BUF_LENGTH) {
68 count = sprintf(buf,
69 "PANIC: printf exceeded %d bytes. Stack has been corrupted.\n",
70 PRINTF_BUF_LENGTH);
71 outwrite(mem, buf, count);
72 }
73 va_end(args);
74 return count;
75 }
76
errprintf(const char * fmt,...)77 int errprintf(const char *fmt, ...)
78 {
79 int count;
80 char buf[PRINTF_BUF_LENGTH];
81 va_list args;
82
83 va_start(args, fmt);
84
85 count = vsprintf(buf, fmt, args);
86 errwrite(buf, count);
87 if (count >= PRINTF_BUF_LENGTH) {
88 count = sprintf(buf,
89 "PANIC: printf exceeded %d bytes. Stack has been corrupted.\n",
90 PRINTF_BUF_LENGTH);
91 errwrite(buf, count);
92 }
93 va_end(args);
94 return count;
95 }
96
97 /* ------ Debugging ------ */
98
99 /* Ghostscript writes debugging output to gs_debug_out. */
100 /* We define gs_debug and gs_debug_out even if DEBUG isn't defined, */
101 /* so that we can compile individual modules with DEBUG set. */
102 char gs_debug[128];
103 FILE *gs_debug_out;
104
105 /* Test whether a given debugging option is selected. */
106 /* Upper-case letters automatically include their lower-case counterpart. */
107 bool
gs_debug_c(int c)108 gs_debug_c(int c)
109 {
110 return
111 (c >= 'a' && c <= 'z' ? gs_debug[c] | gs_debug[c ^ 32] : gs_debug[c]);
112 }
113
114 /* Define the formats for debugging printout. */
115 const char *const dprintf_file_and_line_format = "%10s(%4d): ";
116 const char *const dprintf_file_only_format = "%10s(unkn): ";
117
118 /*
119 * Define the trace printout procedures. We always include these, in case
120 * other modules were compiled with DEBUG set. Note that they must use
121 * out/errprintf, not fprintf nor fput[cs], because of the way that
122 * stdout/stderr are implemented on DLL/shared library builds.
123 */
124 void
dflush(void)125 dflush(void)
126 {
127 errflush();
128 }
129 private const char *
dprintf_file_tail(const char * file)130 dprintf_file_tail(const char *file)
131 {
132 const char *tail = file + strlen(file);
133
134 while (tail > file &&
135 (isalnum(tail[-1]) || tail[-1] == '.' || tail[-1] == '_')
136 )
137 --tail;
138 return tail;
139 }
140 #if __LINE__ /* compiler provides it */
141 void
dprintf_file_and_line(const char * file,int line)142 dprintf_file_and_line(const char *file, int line)
143 {
144 if (gs_debug['/'])
145 dpf(dprintf_file_and_line_format,
146 dprintf_file_tail(file), line);
147 }
148 #else
149 void
dprintf_file_only(const char * file)150 dprintf_file_only(const char *file)
151 {
152 if (gs_debug['/'])
153 dpf(dprintf_file_only_format, dprintf_file_tail(file));
154 }
155 #endif
156 void
printf_program_ident(const gs_memory_t * mem,const char * program_name,long revision_number)157 printf_program_ident(const gs_memory_t *mem, const char *program_name, long revision_number)
158 {
159 if (program_name)
160 outprintf(mem, (revision_number ? "%s " : "%s"), program_name);
161 if (revision_number) {
162 int fpart = revision_number % 100;
163
164 outprintf(mem, "%d.%02d", (int)(revision_number / 100), fpart);
165 }
166 }
167 void
eprintf_program_ident(const char * program_name,long revision_number)168 eprintf_program_ident(const char *program_name,
169 long revision_number)
170 {
171 if (program_name) {
172 epf((revision_number ? "%s " : "%s"), program_name);
173 if (revision_number) {
174 int fpart = revision_number % 100;
175
176 epf("%d.%02d", (int)(revision_number / 100), fpart);
177 }
178 epf(": ");
179 }
180 }
181 #if __LINE__ /* compiler provides it */
182 void
lprintf_file_and_line(const char * file,int line)183 lprintf_file_and_line(const char *file, int line)
184 {
185 epf("%s(%d): ", file, line);
186 }
187 #else
188 void
lprintf_file_only(FILE * f,const char * file)189 lprintf_file_only(FILE * f, const char *file)
190 {
191 epf("%s(?): ", file);
192 }
193 #endif
194
195 /* Log an error return. We always include this, in case other */
196 /* modules were compiled with DEBUG set. */
197 #undef gs_log_error /* in case DEBUG isn't set */
198 int
gs_log_error(int err,const char * file,int line)199 gs_log_error(int err, const char *file, int line)
200 {
201 if (gs_log_errors) {
202 if (file == NULL)
203 dprintf1("Returning error %d.\n", err);
204 else
205 dprintf3("%s(%d): Returning error %d.\n",
206 (const char *)file, line, err);
207 }
208 return err;
209 }
210
211 /* Check for interrupts before a return. */
212 int
gs_return_check_interrupt(const gs_memory_t * mem,int code)213 gs_return_check_interrupt(const gs_memory_t *mem, int code)
214 {
215 if (code < 0)
216 return code;
217 {
218 int icode = gp_check_interrupts(mem);
219
220 return (icode == 0 ? code :
221 gs_note_error((icode > 0 ? gs_error_interrupt : icode)));
222 }
223 }
224
225 /* ------ Substitutes for missing C library functions ------ */
226
227 #ifdef MEMORY__NEED_MEMMOVE /* see memory_.h */
228 /* Copy bytes like memcpy, guaranteed to handle overlap correctly. */
229 /* ANSI C defines the returned value as being the src argument, */
230 /* but with the const restriction removed! */
231 void *
gs_memmove(void * dest,const void * src,size_t len)232 gs_memmove(void *dest, const void *src, size_t len)
233 {
234 if (!len)
235 return (void *)src;
236 #define bdest ((byte *)dest)
237 #define bsrc ((const byte *)src)
238 /* We use len-1 for comparisons because adding len */
239 /* might produce an offset overflow on segmented systems. */
240 if (PTR_LE(bdest, bsrc)) {
241 register byte *end = bdest + (len - 1);
242
243 if (PTR_LE(bsrc, end)) {
244 /* Source overlaps destination from above. */
245 register const byte *from = bsrc;
246 register byte *to = bdest;
247
248 for (;;) {
249 *to = *from;
250 if (to >= end) /* faster than = */
251 return (void *)src;
252 to++;
253 from++;
254 }
255 }
256 } else {
257 register const byte *from = bsrc + (len - 1);
258
259 if (PTR_LE(bdest, from)) {
260 /* Source overlaps destination from below. */
261 register const byte *end = bsrc;
262 register byte *to = bdest + (len - 1);
263
264 for (;;) {
265 *to = *from;
266 if (from <= end) /* faster than = */
267 return (void *)src;
268 to--;
269 from--;
270 }
271 }
272 }
273 #undef bdest
274 #undef bsrc
275 /* No overlap, it's safe to use memcpy. */
276 memcpy(dest, src, len);
277 return (void *)src;
278 }
279 #endif
280
281 #ifdef MEMORY__NEED_MEMCPY /* see memory_.h */
282 void *
gs_memcpy(void * dest,const void * src,size_t len)283 gs_memcpy(void *dest, const void *src, size_t len)
284 {
285 if (len > 0) {
286 #define bdest ((byte *)dest)
287 #define bsrc ((const byte *)src)
288 /* We can optimize this much better later on. */
289 register byte *end = bdest + (len - 1);
290 register const byte *from = bsrc;
291 register byte *to = bdest;
292
293 for (;;) {
294 *to = *from;
295 if (to >= end) /* faster than = */
296 break;
297 to++;
298 from++;
299 }
300 }
301 #undef bdest
302 #undef bsrc
303 return (void *)src;
304 }
305 #endif
306
307 #ifdef MEMORY__NEED_MEMCHR /* see memory_.h */
308 /* ch should obviously be char rather than int, */
309 /* but the ANSI standard declaration uses int. */
310 void *
gs_memchr(const void * ptr,int ch,size_t len)311 gs_memchr(const void *ptr, int ch, size_t len)
312 {
313 if (len > 0) {
314 register const char *p = ptr;
315 register uint count = len;
316
317 do {
318 if (*p == (char)ch)
319 return (void *)p;
320 p++;
321 } while (--count);
322 }
323 return 0;
324 }
325 #endif
326
327 #ifdef MEMORY__NEED_MEMSET /* see memory_.h */
328 /* ch should obviously be char rather than int, */
329 /* but the ANSI standard declaration uses int. */
330 void *
gs_memset(void * dest,register int ch,size_t len)331 gs_memset(void *dest, register int ch, size_t len)
332 {
333 /*
334 * This procedure is used a lot to fill large regions of images,
335 * so we take some trouble to optimize it.
336 */
337 register char *p = dest;
338 register size_t count = len;
339
340 ch &= 255;
341 if (len >= sizeof(long) * 3) {
342 long wd = (ch << 24) | (ch << 16) | (ch << 8) | ch;
343
344 while (ALIGNMENT_MOD(p, sizeof(long)))
345 *p++ = (char)ch, --count;
346 for (; count >= sizeof(long) * 4;
347 p += sizeof(long) * 4, count -= sizeof(long) * 4
348 )
349 ((long *)p)[3] = ((long *)p)[2] = ((long *)p)[1] =
350 ((long *)p)[0] = wd;
351 switch (count >> ARCH_LOG2_SIZEOF_LONG) {
352 case 3:
353 *((long *)p) = wd; p += sizeof(long);
354 case 2:
355 *((long *)p) = wd; p += sizeof(long);
356 case 1:
357 *((long *)p) = wd; p += sizeof(long);
358 count &= sizeof(long) - 1;
359 case 0:
360 default: /* can't happen */
361 DO_NOTHING;
362 }
363 }
364 /* Do any leftover bytes. */
365 for (; count > 0; --count)
366 *p++ = (char)ch;
367 return dest;
368 }
369 #endif
370
371 #ifdef malloc__need_realloc /* see malloc_.h */
372 /* Some systems have non-working implementations of realloc. */
373 void *
gs_realloc(void * old_ptr,size_t old_size,size_t new_size)374 gs_realloc(void *old_ptr, size_t old_size, size_t new_size)
375 {
376 void *new_ptr;
377
378 if (new_size) {
379 new_ptr = malloc(new_size);
380 if (new_ptr == NULL)
381 return NULL;
382 } else
383 new_ptr = NULL;
384 /* We have to pass in the old size, since we have no way to */
385 /* determine it otherwise. */
386 if (old_ptr != NULL) {
387 if (new_ptr != NULL)
388 memcpy(new_ptr, old_ptr, min(old_size, new_size));
389 free(old_ptr);
390 }
391 return new_ptr;
392 }
393 #endif
394
395 /* ------ Debugging support ------ */
396
397 /* Dump a region of memory. */
398 void
debug_dump_bytes(const byte * from,const byte * to,const char * msg)399 debug_dump_bytes(const byte * from, const byte * to, const char *msg)
400 {
401 const byte *p = from;
402
403 if (from < to && msg)
404 dprintf1("%s:\n", msg);
405 while (p != to) {
406 const byte *q = min(p + 16, to);
407
408 dprintf1("0x%lx:", (ulong) p);
409 while (p != q)
410 dprintf1(" %02x", *p++);
411 dputc('\n');
412 }
413 }
414
415 /* Dump a bitmap. */
416 void
debug_dump_bitmap(const byte * bits,uint raster,uint height,const char * msg)417 debug_dump_bitmap(const byte * bits, uint raster, uint height, const char *msg)
418 {
419 uint y;
420 const byte *data = bits;
421
422 for (y = 0; y < height; ++y, data += raster)
423 debug_dump_bytes(data, data + raster, (y == 0 ? msg : NULL));
424 }
425
426 /* Print a string. */
427 void
debug_print_string(const byte * chrs,uint len)428 debug_print_string(const byte * chrs, uint len)
429 {
430 uint i;
431
432 for (i = 0; i < len; i++)
433 dputc(chrs[i]);
434 dflush();
435 }
436
437 /* Print a string in hexdump format. */
438 void
debug_print_string_hex(const byte * chrs,uint len)439 debug_print_string_hex(const byte * chrs, uint len)
440 {
441 uint i;
442
443 for (i = 0; i < len; i++)
444 dprintf1("%02x", chrs[i]);
445 dflush();
446 }
447
448 /*
449 * The following code prints a hex stack backtrace on Linux/Intel systems.
450 * It is here to be patched into places where we need to print such a trace
451 * because of gdb's inability to put breakpoints in dynamically created
452 * threads.
453 *
454 * first_arg is the first argument of the procedure into which this code
455 * is patched.
456 */
457 #define BACKTRACE(first_arg)\
458 BEGIN\
459 ulong *fp_ = (ulong *)&first_arg - 2;\
460 for (; fp_ && (fp_[1] & 0xff000000) == 0x08000000; fp_ = (ulong *)*fp_)\
461 dprintf2(" fp=0x%lx ip=0x%lx\n", (ulong)fp_, fp_[1]);\
462 END
463
464 /* ------ Arithmetic ------ */
465
466 /* Compute M modulo N. Requires N > 0; guarantees 0 <= imod(M,N) < N, */
467 /* regardless of the whims of the % operator for negative operands. */
468 int
imod(int m,int n)469 imod(int m, int n)
470 {
471 if (n <= 0)
472 return 0; /* sanity check */
473 if (m >= 0)
474 return m % n;
475 {
476 int r = -m % n;
477
478 return (r == 0 ? 0 : n - r);
479 }
480 }
481
482 /* Compute the GCD of two integers. */
483 int
igcd(int x,int y)484 igcd(int x, int y)
485 {
486 int c = x, d = y;
487
488 if (c < 0)
489 c = -c;
490 if (d < 0)
491 d = -d;
492 while (c != 0 && d != 0)
493 if (c > d)
494 c %= d;
495 else
496 d %= c;
497 return d + c; /* at most one is non-zero */
498 }
499
500 /* Compute X such that A*X = B mod M. See gxarith.h for details. */
501 int
idivmod(int a,int b,int m)502 idivmod(int a, int b, int m)
503 {
504 /*
505 * Use the approach indicated in Knuth vol. 2, section 4.5.2, Algorithm
506 * X (p. 302) and exercise 15 (p. 315, solution p. 523).
507 */
508 int u1 = 0, u3 = m;
509 int v1 = 1, v3 = a;
510 /*
511 * The following loop will terminate with a * u1 = gcd(a, m) mod m.
512 * Then x = u1 * b / gcd(a, m) mod m. Since we require that
513 * gcd(a, m) | gcd(a, b), it follows that gcd(a, m) | b, so the
514 * division is exact.
515 */
516 while (v3) {
517 int q = u3 / v3, t;
518
519 t = u1 - v1 * q, u1 = v1, v1 = t;
520 t = u3 - v3 * q, u3 = v3, v3 = t;
521 }
522 return imod(u1 * b / igcd(a, m), m);
523 }
524
525 /* Compute floor(log2(N)). Requires N > 0. */
526 int
ilog2(int n)527 ilog2(int n)
528 {
529 int m = n, l = 0;
530
531 while (m >= 16)
532 m >>= 4, l += 4;
533 return
534 (m <= 1 ? l :
535 "\000\000\001\001\002\002\002\002\003\003\003\003\003\003\003\003"[m] + l);
536 }
537
538 #if defined(NEED_SET_FMUL2FIXED) && !USE_ASM
539
540 /*
541 * Floating multiply with fixed result, for avoiding floating point in
542 * common coordinate transformations. Assumes IEEE representation,
543 * 16-bit short, 32-bit long. Optimized for the case where the first
544 * operand has no more than 16 mantissa bits, e.g., where it is a user space
545 * coordinate (which are often integers).
546 *
547 * The assembly language version of this code is actually faster than
548 * the FPU, if the code is compiled with FPU_TYPE=0 (which requires taking
549 * a trap on every FPU operation). If there is no FPU, the assembly
550 * language version of this code is over 10 times as fast as the emulated FPU.
551 */
552 int
set_fmul2fixed_(fixed * pr,long a,long b)553 set_fmul2fixed_(fixed * pr, long /*float */ a, long /*float */ b)
554 {
555 ulong ma = (ushort)(a >> 8) | 0x8000;
556 ulong mb = (ushort)(b >> 8) | 0x8000;
557 int e = 260 + _fixed_shift - ( ((byte)(a >> 23)) + ((byte)(b >> 23)) );
558 ulong p1 = ma * (b & 0xff);
559 ulong p = ma * mb;
560
561 #define p_bits (size_of(p) * 8)
562
563 if ((byte) a) { /* >16 mantissa bits */
564 ulong p2 = (a & 0xff) * mb;
565
566 p += ((((uint) (byte) a * (uint) (byte) b) >> 8) + p1 + p2) >> 8;
567 } else
568 p += p1 >> 8;
569 if ((uint) e < p_bits) /* e = -1 is possible */
570 p >>= e;
571 else if (e >= p_bits) { /* also detects a=0 or b=0 */
572 *pr = fixed_0;
573 return 0;
574 } else if (e >= -(p_bits - 1) || p >= 1L << (p_bits - 1 + e))
575 return_error(gs_error_limitcheck);
576 else
577 p <<= -e;
578 *pr = ((a ^ b) < 0 ? -p : p);
579 return 0;
580 }
581 int
set_dfmul2fixed_(fixed * pr,ulong xalo,long b,long xahi)582 set_dfmul2fixed_(fixed * pr, ulong /*double lo */ xalo, long /*float */ b, long /*double hi */ xahi)
583 {
584 return set_fmul2fixed_(pr,
585 (xahi & (3L << 30)) +
586 ((xahi << 3) & 0x3ffffff8) +
587 (xalo >> 29),
588 b);
589 }
590
591 #endif /* NEED_SET_FMUL2FIXED */
592
593 #if USE_FPU_FIXED
594
595 /*
596 * Convert from floating point to fixed point with scaling.
597 * These are efficient algorithms for FPU-less machines.
598 */
599 #define mbits_float 23
600 #define mbits_double 20
601 int
set_float2fixed_(fixed * pr,long vf,int frac_bits)602 set_float2fixed_(fixed * pr, long /*float */ vf, int frac_bits)
603 {
604 fixed mantissa;
605 int shift;
606
607 if (!(vf & 0x7f800000)) {
608 *pr = fixed_0;
609 return 0;
610 }
611 mantissa = (fixed) ((vf & 0x7fffff) | 0x800000);
612 shift = ((vf >> 23) & 255) - (127 + 23) + frac_bits;
613 if (shift >= 0) {
614 if (shift >= sizeof(fixed) * 8 - 24)
615 return_error(gs_error_limitcheck);
616 if (vf < 0)
617 mantissa = -mantissa;
618 *pr = (fixed) (mantissa << shift);
619 } else
620 *pr = (shift < -24 ? fixed_0 :
621 vf < 0 ? -(fixed) (mantissa >> -shift) : /* truncate */
622 (fixed) (mantissa >> -shift));
623 return 0;
624 }
625 int
set_double2fixed_(fixed * pr,ulong lo,long hi,int frac_bits)626 set_double2fixed_(fixed * pr, ulong /*double lo */ lo,
627 long /*double hi */ hi, int frac_bits)
628 {
629 fixed mantissa;
630 int shift;
631
632 if (!(hi & 0x7ff00000)) {
633 *pr = fixed_0;
634 return 0;
635 }
636 /* We only use 31 bits of mantissa even if sizeof(long) > 4. */
637 mantissa = (fixed) (((hi & 0xfffff) << 10) | (lo >> 22) | 0x40000000);
638 shift = ((hi >> 20) & 2047) - (1023 + 30) + frac_bits;
639 if (shift > 0)
640 return_error(gs_error_limitcheck);
641 *pr = (shift < -30 ? fixed_0 :
642 hi < 0 ? -(fixed) (mantissa >> -shift) : /* truncate */
643 (fixed) (mantissa >> -shift));
644 return 0;
645 }
646 /*
647 * Given a fixed value x with fbits bits of fraction, set v to the mantissa
648 * (left-justified in 32 bits) and f to the exponent word of the
649 * corresponding floating-point value with mbits bits of mantissa in the
650 * first word. (The exponent part of f is biased by -1, because we add the
651 * top 1-bit of the mantissa to it.)
652 */
653 static const byte f2f_shifts[] =
654 {4, 3, 2, 2, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0};
655
656 #define f2f_declare(v, f)\
657 ulong v;\
658 long f
659 #define f2f(x, v, f, mbits, fbits)\
660 if ( x < 0 )\
661 f = 0xc0000000 + (29 << mbits) - ((long)fbits << mbits), v = -x;\
662 else\
663 f = 0x40000000 + (29 << mbits) - ((long)fbits << mbits), v = x;\
664 if ( v < 0x8000 )\
665 v <<= 15, f -= 15 << mbits;\
666 if ( v < 0x800000 )\
667 v <<= 8, f -= 8 << mbits;\
668 if ( v < 0x8000000 )\
669 v <<= 4, f -= 4 << mbits;\
670 { int shift = f2f_shifts[v >> 28];\
671 v <<= shift, f -= shift << mbits;\
672 }
673 long
fixed2float_(fixed x,int frac_bits)674 fixed2float_(fixed x, int frac_bits)
675 {
676 f2f_declare(v, f);
677
678 if (x == 0)
679 return 0;
680 f2f(x, v, f, mbits_float, frac_bits);
681 return f + (((v >> 7) + 1) >> 1);
682 }
683 void
set_fixed2double_(double * pd,fixed x,int frac_bits)684 set_fixed2double_(double *pd, fixed x, int frac_bits)
685 {
686 f2f_declare(v, f);
687
688 if (x == 0) {
689 ((long *)pd)[1 - arch_is_big_endian] = 0;
690 ((ulong *) pd)[arch_is_big_endian] = 0;
691 } else {
692 f2f(x, v, f, mbits_double, frac_bits);
693 ((long *)pd)[1 - arch_is_big_endian] = f + (v >> 11);
694 ((ulong *) pd)[arch_is_big_endian] = v << 21;
695 }
696 }
697
698 #endif /* USE_FPU_FIXED */
699
700 /*
701 * Compute A * B / C when 0 <= B < C and A * B exceeds (or might exceed)
702 * the capacity of a long.
703 * Note that this procedure takes the floor, rather than truncating
704 * towards zero, if A < 0. This ensures that 0 <= R < C.
705 */
706
707 #define num_bits (sizeof(fixed) * 8)
708 #define half_bits (num_bits / 2)
709 #define half_mask ((1L << half_bits) - 1)
710
711 /*
712 * If doubles aren't wide enough, we lose too much precision by using double
713 * arithmetic: we have to use the slower, accurate fixed-point algorithm.
714 * See the simpler implementation below for more information.
715 */
716 #define MAX_OTHER_FACTOR_BITS\
717 (ARCH_DOUBLE_MANTISSA_BITS - ARCH_SIZEOF_FIXED * 8)
718 #define ROUND_BITS\
719 (ARCH_SIZEOF_FIXED * 8 * 2 - ARCH_DOUBLE_MANTISSA_BITS)
720
721 #if USE_FPU_FIXED || ROUND_BITS >= MAX_OTHER_FACTOR_BITS - 1
722
723 #ifdef DEBUG
724 struct {
725 long mnanb, mnab, manb, mab, mnc, mdq, mde, mds, mqh, mql;
726 } fmq_stat;
727 # define mincr(x) ++fmq_stat.x
728 #else
729 # define mincr(x) DO_NOTHING
730 #endif
731 fixed
fixed_mult_quo(fixed signed_A,fixed B,fixed C)732 fixed_mult_quo(fixed signed_A, fixed B, fixed C)
733 {
734 /* First compute A * B in double-fixed precision. */
735 ulong A = (signed_A < 0 ? -signed_A : signed_A);
736 long msw;
737 ulong lsw;
738 ulong p1;
739
740 if (B <= half_mask) {
741 if (A <= half_mask) {
742 ulong P = A * B;
743 fixed Q = P / (ulong)C;
744
745 mincr(mnanb);
746 /* If A < 0 and the division isn't exact, take the floor. */
747 return (signed_A >= 0 ? Q : Q * C == P ? -Q : ~Q /* -Q - 1 */);
748 }
749 /*
750 * We might still have C <= half_mask, which we can
751 * handle with a simpler algorithm.
752 */
753 lsw = (A & half_mask) * B;
754 p1 = (A >> half_bits) * B;
755 if (C <= half_mask) {
756 fixed q0 = (p1 += lsw >> half_bits) / C;
757 ulong rem = ((p1 - C * q0) << half_bits) + (lsw & half_mask);
758 ulong q1 = rem / (ulong)C;
759 fixed Q = (q0 << half_bits) + q1;
760
761 mincr(mnc);
762 /* If A < 0 and the division isn't exact, take the floor. */
763 return (signed_A >= 0 ? Q : q1 * C == rem ? -Q : ~Q);
764 }
765 msw = p1 >> half_bits;
766 mincr(manb);
767 } else if (A <= half_mask) {
768 p1 = A * (B >> half_bits);
769 msw = p1 >> half_bits;
770 lsw = A * (B & half_mask);
771 mincr(mnab);
772 } else { /* We have to compute all 4 products. :-( */
773 ulong lo_A = A & half_mask;
774 ulong hi_A = A >> half_bits;
775 ulong lo_B = B & half_mask;
776 ulong hi_B = B >> half_bits;
777 ulong p1x = hi_A * lo_B;
778
779 msw = hi_A * hi_B;
780 lsw = lo_A * lo_B;
781 p1 = lo_A * hi_B;
782 if (p1 > max_ulong - p1x)
783 msw += 1L << half_bits;
784 p1 += p1x;
785 msw += p1 >> half_bits;
786 mincr(mab);
787 }
788 /* Finish up by adding the low half of p1 to the high half of lsw. */
789 #if max_fixed < max_long
790 p1 &= half_mask;
791 #endif
792 p1 <<= half_bits;
793 if (p1 > max_ulong - lsw)
794 msw++;
795 lsw += p1;
796 /*
797 * Now divide the double-length product by C. Note that we know msw
798 * < C (otherwise the quotient would overflow). Start by shifting
799 * (msw,lsw) and C left until C >= 1 << (num_bits - 1).
800 */
801 {
802 ulong denom = C;
803 int shift = 0;
804
805 #define bits_4th (num_bits / 4)
806 if (denom < 1L << (num_bits - bits_4th)) {
807 mincr(mdq);
808 denom <<= bits_4th, shift += bits_4th;
809 }
810 #undef bits_4th
811 #define bits_8th (num_bits / 8)
812 if (denom < 1L << (num_bits - bits_8th)) {
813 mincr(mde);
814 denom <<= bits_8th, shift += bits_8th;
815 }
816 #undef bits_8th
817 while (!(denom & (-1L << (num_bits - 1)))) {
818 mincr(mds);
819 denom <<= 1, ++shift;
820 }
821 msw = (msw << shift) + (lsw >> (num_bits - shift));
822 lsw <<= shift;
823 #if max_fixed < max_long
824 lsw &= (1L << (sizeof(fixed) * 8)) - 1;
825 #endif
826 /* Compute a trial upper-half quotient. */
827 {
828 ulong hi_D = denom >> half_bits;
829 ulong lo_D = denom & half_mask;
830 ulong hi_Q = (ulong) msw / hi_D;
831
832 /* hi_Q might be too high by 1 or 2, but it isn't too low. */
833 ulong p0 = hi_Q * hi_D;
834 ulong p1 = hi_Q * lo_D;
835 ulong hi_P;
836
837 while ((hi_P = p0 + (p1 >> half_bits)) > msw ||
838 (hi_P == msw && ((p1 & half_mask) << half_bits) > lsw)
839 ) { /* hi_Q was too high by 1. */
840 --hi_Q;
841 p0 -= hi_D;
842 p1 -= lo_D;
843 mincr(mqh);
844 }
845 p1 = (p1 & half_mask) << half_bits;
846 if (p1 > lsw)
847 msw--;
848 lsw -= p1;
849 msw -= hi_P;
850 /* Now repeat to get the lower-half quotient. */
851 msw = (msw << half_bits) + (lsw >> half_bits);
852 #if max_fixed < max_long
853 lsw &= half_mask;
854 #endif
855 lsw <<= half_bits;
856 {
857 ulong lo_Q = (ulong) msw / hi_D;
858 long Q;
859
860 p1 = lo_Q * lo_D;
861 p0 = lo_Q * hi_D;
862 while ((hi_P = p0 + (p1 >> half_bits)) > msw ||
863 (hi_P == msw && ((p1 & half_mask) << half_bits) > lsw)
864 ) { /* lo_Q was too high by 1. */
865 --lo_Q;
866 p0 -= hi_D;
867 p1 -= lo_D;
868 mincr(mql);
869 }
870 Q = (hi_Q << half_bits) + lo_Q;
871 return (signed_A >= 0 ? Q : p0 | p1 ? ~Q /* -Q - 1 */ : -Q);
872 }
873 }
874 }
875 }
876
877 #else /* use doubles */
878
879 /*
880 * Compute A * B / C as above using doubles. If floating point is
881 * reasonably fast, this is much faster than the fixed-point algorithm.
882 */
883 fixed
fixed_mult_quo(fixed signed_A,fixed B,fixed C)884 fixed_mult_quo(fixed signed_A, fixed B, fixed C)
885 {
886 /*
887 * Check whether A * B will fit in the mantissa of a double.
888 */
889 #define MAX_OTHER_FACTOR (1L << MAX_OTHER_FACTOR_BITS)
890 if (B < MAX_OTHER_FACTOR || any_abs(signed_A) < MAX_OTHER_FACTOR) {
891 #undef MAX_OTHER_FACTOR
892 /*
893 * The product fits, so a straightforward double computation
894 * will be exact.
895 */
896 return (fixed)floor((double)signed_A * B / C);
897 } else {
898 /*
899 * The product won't fit. However, the approximate product will
900 * only be off by at most +/- 1/2 * (1 << ROUND_BITS) because of
901 * rounding. If we add 1 << ROUND_BITS to the value of the product
902 * (i.e., 1 in the least significant bit of the mantissa), the
903 * result is always greater than the correct product by between 1/2
904 * and 3/2 * (1 << ROUND_BITS). We know this is less than C:
905 * because of the 'if' just above, we know that B >=
906 * MAX_OTHER_FACTOR; since B <= C, we know C >= MAX_OTHER_FACTOR;
907 * and because of the #if that chose between the two
908 * implementations, we know that C >= 2 * (1 << ROUND_BITS). Hence,
909 * the quotient after dividing by C will be at most 1 too large.
910 */
911 fixed q =
912 (fixed)floor(((double)signed_A * B + (1L << ROUND_BITS)) / C);
913
914 /*
915 * Compute the remainder R. If the quotient was correct,
916 * 0 <= R < C. If the quotient was too high, -C <= R < 0.
917 */
918 if (signed_A * B - q * C < 0)
919 --q;
920 return q;
921 }
922 }
923
924 #endif
925
926 #undef MAX_OTHER_FACTOR_BITS
927 #undef ROUND_BITS
928
929 #undef num_bits
930 #undef half_bits
931 #undef half_mask
932
933 /* Trace calls on sqrt when debugging. */
934 double
gs_sqrt(double x,const char * file,int line)935 gs_sqrt(double x, const char *file, int line)
936 {
937 if (gs_debug_c('~')) {
938 dprintf3("[~]sqrt(%g) at %s:%d\n", x, (const char *)file, line);
939 dflush();
940 }
941 return orig_sqrt(x);
942 }
943
944 /*
945 * Define sine and cosine functions that take angles in degrees rather than
946 * radians, and that are implemented efficiently on machines with slow
947 * (or no) floating point.
948 */
949 #if USE_FPU < 0 /****** maybe should be <= 0 ? ***** */
950
951 #define sin0 0.00000000000000000
952 #define sin1 0.01745240643728351
953 #define sin2 0.03489949670250097
954 #define sin3 0.05233595624294383
955 #define sin4 0.06975647374412530
956 #define sin5 0.08715574274765817
957 #define sin6 0.10452846326765346
958 #define sin7 0.12186934340514748
959 #define sin8 0.13917310096006544
960 #define sin9 0.15643446504023087
961 #define sin10 0.17364817766693033
962 #define sin11 0.19080899537654480
963 #define sin12 0.20791169081775931
964 #define sin13 0.22495105434386498
965 #define sin14 0.24192189559966773
966 #define sin15 0.25881904510252074
967 #define sin16 0.27563735581699916
968 #define sin17 0.29237170472273671
969 #define sin18 0.30901699437494740
970 #define sin19 0.32556815445715670
971 #define sin20 0.34202014332566871
972 #define sin21 0.35836794954530027
973 #define sin22 0.37460659341591201
974 #define sin23 0.39073112848927377
975 #define sin24 0.40673664307580015
976 #define sin25 0.42261826174069944
977 #define sin26 0.43837114678907740
978 #define sin27 0.45399049973954675
979 #define sin28 0.46947156278589081
980 #define sin29 0.48480962024633706
981 #define sin30 0.50000000000000000
982 #define sin31 0.51503807491005416
983 #define sin32 0.52991926423320490
984 #define sin33 0.54463903501502708
985 #define sin34 0.55919290347074679
986 #define sin35 0.57357643635104605
987 #define sin36 0.58778525229247314
988 #define sin37 0.60181502315204827
989 #define sin38 0.61566147532565829
990 #define sin39 0.62932039104983739
991 #define sin40 0.64278760968653925
992 #define sin41 0.65605902899050728
993 #define sin42 0.66913060635885824
994 #define sin43 0.68199836006249848
995 #define sin44 0.69465837045899725
996 #define sin45 0.70710678118654746
997 #define sin46 0.71933980033865108
998 #define sin47 0.73135370161917046
999 #define sin48 0.74314482547739413
1000 #define sin49 0.75470958022277201
1001 #define sin50 0.76604444311897801
1002 #define sin51 0.77714596145697090
1003 #define sin52 0.78801075360672190
1004 #define sin53 0.79863551004729283
1005 #define sin54 0.80901699437494745
1006 #define sin55 0.81915204428899180
1007 #define sin56 0.82903757255504174
1008 #define sin57 0.83867056794542394
1009 #define sin58 0.84804809615642596
1010 #define sin59 0.85716730070211222
1011 #define sin60 0.86602540378443860
1012 #define sin61 0.87461970713939574
1013 #define sin62 0.88294759285892688
1014 #define sin63 0.89100652418836779
1015 #define sin64 0.89879404629916704
1016 #define sin65 0.90630778703664994
1017 #define sin66 0.91354545764260087
1018 #define sin67 0.92050485345244037
1019 #define sin68 0.92718385456678731
1020 #define sin69 0.93358042649720174
1021 #define sin70 0.93969262078590832
1022 #define sin71 0.94551857559931674
1023 #define sin72 0.95105651629515353
1024 #define sin73 0.95630475596303544
1025 #define sin74 0.96126169593831889
1026 #define sin75 0.96592582628906831
1027 #define sin76 0.97029572627599647
1028 #define sin77 0.97437006478523525
1029 #define sin78 0.97814760073380558
1030 #define sin79 0.98162718344766398
1031 #define sin80 0.98480775301220802
1032 #define sin81 0.98768834059513777
1033 #define sin82 0.99026806874157036
1034 #define sin83 0.99254615164132198
1035 #define sin84 0.99452189536827329
1036 #define sin85 0.99619469809174555
1037 #define sin86 0.99756405025982420
1038 #define sin87 0.99862953475457383
1039 #define sin88 0.99939082701909576
1040 #define sin89 0.99984769515639127
1041 #define sin90 1.00000000000000000
1042
1043 private const double sin_table[361] =
1044 {
1045 sin0,
1046 sin1, sin2, sin3, sin4, sin5, sin6, sin7, sin8, sin9, sin10,
1047 sin11, sin12, sin13, sin14, sin15, sin16, sin17, sin18, sin19, sin20,
1048 sin21, sin22, sin23, sin24, sin25, sin26, sin27, sin28, sin29, sin30,
1049 sin31, sin32, sin33, sin34, sin35, sin36, sin37, sin38, sin39, sin40,
1050 sin41, sin42, sin43, sin44, sin45, sin46, sin47, sin48, sin49, sin50,
1051 sin51, sin52, sin53, sin54, sin55, sin56, sin57, sin58, sin59, sin60,
1052 sin61, sin62, sin63, sin64, sin65, sin66, sin67, sin68, sin69, sin70,
1053 sin71, sin72, sin73, sin74, sin75, sin76, sin77, sin78, sin79, sin80,
1054 sin81, sin82, sin83, sin84, sin85, sin86, sin87, sin88, sin89, sin90,
1055 sin89, sin88, sin87, sin86, sin85, sin84, sin83, sin82, sin81, sin80,
1056 sin79, sin78, sin77, sin76, sin75, sin74, sin73, sin72, sin71, sin70,
1057 sin69, sin68, sin67, sin66, sin65, sin64, sin63, sin62, sin61, sin60,
1058 sin59, sin58, sin57, sin56, sin55, sin54, sin53, sin52, sin51, sin50,
1059 sin49, sin48, sin47, sin46, sin45, sin44, sin43, sin42, sin41, sin40,
1060 sin39, sin38, sin37, sin36, sin35, sin34, sin33, sin32, sin31, sin30,
1061 sin29, sin28, sin27, sin26, sin25, sin24, sin23, sin22, sin21, sin20,
1062 sin19, sin18, sin17, sin16, sin15, sin14, sin13, sin12, sin11, sin10,
1063 sin9, sin8, sin7, sin6, sin5, sin4, sin3, sin2, sin1, sin0,
1064 -sin1, -sin2, -sin3, -sin4, -sin5, -sin6, -sin7, -sin8, -sin9, -sin10,
1065 -sin11, -sin12, -sin13, -sin14, -sin15, -sin16, -sin17, -sin18, -sin19, -sin20,
1066 -sin21, -sin22, -sin23, -sin24, -sin25, -sin26, -sin27, -sin28, -sin29, -sin30,
1067 -sin31, -sin32, -sin33, -sin34, -sin35, -sin36, -sin37, -sin38, -sin39, -sin40,
1068 -sin41, -sin42, -sin43, -sin44, -sin45, -sin46, -sin47, -sin48, -sin49, -sin50,
1069 -sin51, -sin52, -sin53, -sin54, -sin55, -sin56, -sin57, -sin58, -sin59, -sin60,
1070 -sin61, -sin62, -sin63, -sin64, -sin65, -sin66, -sin67, -sin68, -sin69, -sin70,
1071 -sin71, -sin72, -sin73, -sin74, -sin75, -sin76, -sin77, -sin78, -sin79, -sin80,
1072 -sin81, -sin82, -sin83, -sin84, -sin85, -sin86, -sin87, -sin88, -sin89, -sin90,
1073 -sin89, -sin88, -sin87, -sin86, -sin85, -sin84, -sin83, -sin82, -sin81, -sin80,
1074 -sin79, -sin78, -sin77, -sin76, -sin75, -sin74, -sin73, -sin72, -sin71, -sin70,
1075 -sin69, -sin68, -sin67, -sin66, -sin65, -sin64, -sin63, -sin62, -sin61, -sin60,
1076 -sin59, -sin58, -sin57, -sin56, -sin55, -sin54, -sin53, -sin52, -sin51, -sin50,
1077 -sin49, -sin48, -sin47, -sin46, -sin45, -sin44, -sin43, -sin42, -sin41, -sin40,
1078 -sin39, -sin38, -sin37, -sin36, -sin35, -sin34, -sin33, -sin32, -sin31, -sin30,
1079 -sin29, -sin28, -sin27, -sin26, -sin25, -sin24, -sin23, -sin22, -sin21, -sin20,
1080 -sin19, -sin18, -sin17, -sin16, -sin15, -sin14, -sin13, -sin12, -sin11, -sin10,
1081 -sin9, -sin8, -sin7, -sin6, -sin5, -sin4, -sin3, -sin2, -sin1, -sin0
1082 };
1083
1084 double
gs_sin_degrees(double ang)1085 gs_sin_degrees(double ang)
1086 {
1087 int ipart;
1088
1089 if (is_fneg(ang))
1090 ang = 180 - ang;
1091 ipart = (int)ang;
1092 if (ipart >= 360) {
1093 int arem = ipart % 360;
1094
1095 ang -= (ipart - arem);
1096 ipart = arem;
1097 }
1098 return
1099 (ang == ipart ? sin_table[ipart] :
1100 sin_table[ipart] + (sin_table[ipart + 1] - sin_table[ipart]) *
1101 (ang - ipart));
1102 }
1103
1104 double
gs_cos_degrees(double ang)1105 gs_cos_degrees(double ang)
1106 {
1107 int ipart;
1108
1109 if (is_fneg(ang))
1110 ang = 90 - ang;
1111 else
1112 ang += 90;
1113 ipart = (int)ang;
1114 if (ipart >= 360) {
1115 int arem = ipart % 360;
1116
1117 ang -= (ipart - arem);
1118 ipart = arem;
1119 }
1120 return
1121 (ang == ipart ? sin_table[ipart] :
1122 sin_table[ipart] + (sin_table[ipart + 1] - sin_table[ipart]) *
1123 (ang - ipart));
1124 }
1125
1126 void
gs_sincos_degrees(double ang,gs_sincos_t * psincos)1127 gs_sincos_degrees(double ang, gs_sincos_t * psincos)
1128 {
1129 psincos->sin = gs_sin_degrees(ang);
1130 psincos->cos = gs_cos_degrees(ang);
1131 psincos->orthogonal =
1132 (is_fzero(psincos->sin) || is_fzero(psincos->cos));
1133 }
1134
1135 #else /* we have floating point */
1136
1137 static const int isincos[5] =
1138 {0, 1, 0, -1, 0};
1139
1140 /* GCC with -ffast-math compiles ang/90. as ang*(1/90.), losing precission.
1141 * This doesn't happen when the numeral is replaced with a non-const variable.
1142 * So we define the variable to work around the GCC problem.
1143 */
1144 static double const_90_degrees = 90.;
1145
1146 double
gs_sin_degrees(double ang)1147 gs_sin_degrees(double ang)
1148 {
1149 double quot = ang / const_90_degrees;
1150
1151 if (floor(quot) == quot) {
1152 /*
1153 * We need 4.0, rather than 4, here because of non-ANSI compilers.
1154 * The & 3 is because quot might be negative.
1155 */
1156 return isincos[(int)fmod(quot, 4.0) & 3];
1157 }
1158 return sin(ang * (M_PI / 180));
1159 }
1160
1161 double
gs_cos_degrees(double ang)1162 gs_cos_degrees(double ang)
1163 {
1164 double quot = ang / const_90_degrees;
1165
1166 if (floor(quot) == quot) {
1167 /* See above re the following line. */
1168 return isincos[((int)fmod(quot, 4.0) & 3) + 1];
1169 }
1170 return cos(ang * (M_PI / 180));
1171 }
1172
1173 void
gs_sincos_degrees(double ang,gs_sincos_t * psincos)1174 gs_sincos_degrees(double ang, gs_sincos_t * psincos)
1175 {
1176 double quot = ang / const_90_degrees;
1177
1178 if (floor(quot) == quot) {
1179 /* See above re the following line. */
1180 int quads = (int)fmod(quot, 4.0) & 3;
1181
1182 psincos->sin = isincos[quads];
1183 psincos->cos = isincos[quads + 1];
1184 psincos->orthogonal = true;
1185 } else {
1186 double arad = ang * (M_PI / 180);
1187
1188 psincos->sin = sin(arad);
1189 psincos->cos = cos(arad);
1190 psincos->orthogonal = false;
1191 }
1192 }
1193
1194 #endif /* USE_FPU */
1195
1196 /*
1197 * Define an atan2 function that returns an angle in degrees and uses
1198 * the PostScript quadrant rules. Note that it may return
1199 * gs_error_undefinedresult.
1200 */
1201 int
gs_atan2_degrees(double y,double x,double * pangle)1202 gs_atan2_degrees(double y, double x, double *pangle)
1203 {
1204 if (y == 0) { /* on X-axis, special case */
1205 if (x == 0)
1206 return_error(gs_error_undefinedresult);
1207 *pangle = (x < 0 ? 180 : 0);
1208 } else {
1209 double result = atan2(y, x) * radians_to_degrees;
1210
1211 if (result < 0)
1212 result += 360;
1213 *pangle = result;
1214 }
1215 return 0;
1216 }
1217