xref: /openbmc/linux/arch/x86/math-emu/fpu_trig.c (revision 6c9111bc)
1 // SPDX-License-Identifier: GPL-2.0
2 /*---------------------------------------------------------------------------+
3  |  fpu_trig.c                                                               |
4  |                                                                           |
5  | Implementation of the FPU "transcendental" functions.                     |
6  |                                                                           |
7  | Copyright (C) 1992,1993,1994,1997,1999                                    |
8  |                       W. Metzenthen, 22 Parker St, Ormond, Vic 3163,      |
9  |                       Australia.  E-mail   billm@melbpc.org.au            |
10  |                                                                           |
11  |                                                                           |
12  +---------------------------------------------------------------------------*/
13 
14 #include "fpu_system.h"
15 #include "exception.h"
16 #include "fpu_emu.h"
17 #include "status_w.h"
18 #include "control_w.h"
19 #include "reg_constant.h"
20 
21 static void rem_kernel(unsigned long long st0, unsigned long long *y,
22 		       unsigned long long st1, unsigned long long q, int n);
23 
24 #define BETTER_THAN_486
25 
26 #define FCOS  4
27 
28 /* Used only by fptan, fsin, fcos, and fsincos. */
29 /* This routine produces very accurate results, similar to
30    using a value of pi with more than 128 bits precision. */
31 /* Limited measurements show no results worse than 64 bit precision
32    except for the results for arguments close to 2^63, where the
33    precision of the result sometimes degrades to about 63.9 bits */
34 static int trig_arg(FPU_REG *st0_ptr, int even)
35 {
36 	FPU_REG tmp;
37 	u_char tmptag;
38 	unsigned long long q;
39 	int old_cw = control_word, saved_status = partial_status;
40 	int tag, st0_tag = TAG_Valid;
41 
42 	if (exponent(st0_ptr) >= 63) {
43 		partial_status |= SW_C2;	/* Reduction incomplete. */
44 		return -1;
45 	}
46 
47 	control_word &= ~CW_RC;
48 	control_word |= RC_CHOP;
49 
50 	setpositive(st0_ptr);
51 	tag = FPU_u_div(st0_ptr, &CONST_PI2, &tmp, PR_64_BITS | RC_CHOP | 0x3f,
52 			SIGN_POS);
53 
54 	FPU_round_to_int(&tmp, tag);	/* Fortunately, this can't overflow
55 					   to 2^64 */
56 	q = significand(&tmp);
57 	if (q) {
58 		rem_kernel(significand(st0_ptr),
59 			   &significand(&tmp),
60 			   significand(&CONST_PI2),
61 			   q, exponent(st0_ptr) - exponent(&CONST_PI2));
62 		setexponent16(&tmp, exponent(&CONST_PI2));
63 		st0_tag = FPU_normalize(&tmp);
64 		FPU_copy_to_reg0(&tmp, st0_tag);
65 	}
66 
67 	if ((even && !(q & 1)) || (!even && (q & 1))) {
68 		st0_tag =
69 		    FPU_sub(REV | LOADED | TAG_Valid, (int)&CONST_PI2,
70 			    FULL_PRECISION);
71 
72 #ifdef BETTER_THAN_486
73 		/* So far, the results are exact but based upon a 64 bit
74 		   precision approximation to pi/2. The technique used
75 		   now is equivalent to using an approximation to pi/2 which
76 		   is accurate to about 128 bits. */
77 		if ((exponent(st0_ptr) <= exponent(&CONST_PI2extra) + 64)
78 		    || (q > 1)) {
79 			/* This code gives the effect of having pi/2 to better than
80 			   128 bits precision. */
81 
82 			significand(&tmp) = q + 1;
83 			setexponent16(&tmp, 63);
84 			FPU_normalize(&tmp);
85 			tmptag =
86 			    FPU_u_mul(&CONST_PI2extra, &tmp, &tmp,
87 				      FULL_PRECISION, SIGN_POS,
88 				      exponent(&CONST_PI2extra) +
89 				      exponent(&tmp));
90 			setsign(&tmp, getsign(&CONST_PI2extra));
91 			st0_tag = FPU_add(&tmp, tmptag, 0, FULL_PRECISION);
92 			if (signnegative(st0_ptr)) {
93 				/* CONST_PI2extra is negative, so the result of the addition
94 				   can be negative. This means that the argument is actually
95 				   in a different quadrant. The correction is always < pi/2,
96 				   so it can't overflow into yet another quadrant. */
97 				setpositive(st0_ptr);
98 				q++;
99 			}
100 		}
101 #endif /* BETTER_THAN_486 */
102 	}
103 #ifdef BETTER_THAN_486
104 	else {
105 		/* So far, the results are exact but based upon a 64 bit
106 		   precision approximation to pi/2. The technique used
107 		   now is equivalent to using an approximation to pi/2 which
108 		   is accurate to about 128 bits. */
109 		if (((q > 0)
110 		     && (exponent(st0_ptr) <= exponent(&CONST_PI2extra) + 64))
111 		    || (q > 1)) {
112 			/* This code gives the effect of having p/2 to better than
113 			   128 bits precision. */
114 
115 			significand(&tmp) = q;
116 			setexponent16(&tmp, 63);
117 			FPU_normalize(&tmp);	/* This must return TAG_Valid */
118 			tmptag =
119 			    FPU_u_mul(&CONST_PI2extra, &tmp, &tmp,
120 				      FULL_PRECISION, SIGN_POS,
121 				      exponent(&CONST_PI2extra) +
122 				      exponent(&tmp));
123 			setsign(&tmp, getsign(&CONST_PI2extra));
124 			st0_tag = FPU_sub(LOADED | (tmptag & 0x0f), (int)&tmp,
125 					  FULL_PRECISION);
126 			if ((exponent(st0_ptr) == exponent(&CONST_PI2)) &&
127 			    ((st0_ptr->sigh > CONST_PI2.sigh)
128 			     || ((st0_ptr->sigh == CONST_PI2.sigh)
129 				 && (st0_ptr->sigl > CONST_PI2.sigl)))) {
130 				/* CONST_PI2extra is negative, so the result of the
131 				   subtraction can be larger than pi/2. This means
132 				   that the argument is actually in a different quadrant.
133 				   The correction is always < pi/2, so it can't overflow
134 				   into yet another quadrant. */
135 				st0_tag =
136 				    FPU_sub(REV | LOADED | TAG_Valid,
137 					    (int)&CONST_PI2, FULL_PRECISION);
138 				q++;
139 			}
140 		}
141 	}
142 #endif /* BETTER_THAN_486 */
143 
144 	FPU_settag0(st0_tag);
145 	control_word = old_cw;
146 	partial_status = saved_status & ~SW_C2;	/* Reduction complete. */
147 
148 	return (q & 3) | even;
149 }
150 
151 /* Convert a long to register */
152 static void convert_l2reg(long const *arg, int deststnr)
153 {
154 	int tag;
155 	long num = *arg;
156 	u_char sign;
157 	FPU_REG *dest = &st(deststnr);
158 
159 	if (num == 0) {
160 		FPU_copy_to_regi(&CONST_Z, TAG_Zero, deststnr);
161 		return;
162 	}
163 
164 	if (num > 0) {
165 		sign = SIGN_POS;
166 	} else {
167 		num = -num;
168 		sign = SIGN_NEG;
169 	}
170 
171 	dest->sigh = num;
172 	dest->sigl = 0;
173 	setexponent16(dest, 31);
174 	tag = FPU_normalize(dest);
175 	FPU_settagi(deststnr, tag);
176 	setsign(dest, sign);
177 	return;
178 }
179 
180 static void single_arg_error(FPU_REG *st0_ptr, u_char st0_tag)
181 {
182 	if (st0_tag == TAG_Empty)
183 		FPU_stack_underflow();	/* Puts a QNaN in st(0) */
184 	else if (st0_tag == TW_NaN)
185 		real_1op_NaN(st0_ptr);	/* return with a NaN in st(0) */
186 #ifdef PARANOID
187 	else
188 		EXCEPTION(EX_INTERNAL | 0x0112);
189 #endif /* PARANOID */
190 }
191 
192 static void single_arg_2_error(FPU_REG *st0_ptr, u_char st0_tag)
193 {
194 	int isNaN;
195 
196 	switch (st0_tag) {
197 	case TW_NaN:
198 		isNaN = (exponent(st0_ptr) == EXP_OVER)
199 		    && (st0_ptr->sigh & 0x80000000);
200 		if (isNaN && !(st0_ptr->sigh & 0x40000000)) {	/* Signaling ? */
201 			EXCEPTION(EX_Invalid);
202 			if (control_word & CW_Invalid) {
203 				/* The masked response */
204 				/* Convert to a QNaN */
205 				st0_ptr->sigh |= 0x40000000;
206 				push();
207 				FPU_copy_to_reg0(st0_ptr, TAG_Special);
208 			}
209 		} else if (isNaN) {
210 			/* A QNaN */
211 			push();
212 			FPU_copy_to_reg0(st0_ptr, TAG_Special);
213 		} else {
214 			/* pseudoNaN or other unsupported */
215 			EXCEPTION(EX_Invalid);
216 			if (control_word & CW_Invalid) {
217 				/* The masked response */
218 				FPU_copy_to_reg0(&CONST_QNaN, TAG_Special);
219 				push();
220 				FPU_copy_to_reg0(&CONST_QNaN, TAG_Special);
221 			}
222 		}
223 		break;		/* return with a NaN in st(0) */
224 #ifdef PARANOID
225 	default:
226 		EXCEPTION(EX_INTERNAL | 0x0112);
227 #endif /* PARANOID */
228 	}
229 }
230 
231 /*---------------------------------------------------------------------------*/
232 
233 static void f2xm1(FPU_REG *st0_ptr, u_char tag)
234 {
235 	FPU_REG a;
236 
237 	clear_C1();
238 
239 	if (tag == TAG_Valid) {
240 		/* For an 80486 FPU, the result is undefined if the arg is >= 1.0 */
241 		if (exponent(st0_ptr) < 0) {
242 		      denormal_arg:
243 
244 			FPU_to_exp16(st0_ptr, &a);
245 
246 			/* poly_2xm1(x) requires 0 < st(0) < 1. */
247 			poly_2xm1(getsign(st0_ptr), &a, st0_ptr);
248 		}
249 		set_precision_flag_up();	/* 80486 appears to always do this */
250 		return;
251 	}
252 
253 	if (tag == TAG_Zero)
254 		return;
255 
256 	if (tag == TAG_Special)
257 		tag = FPU_Special(st0_ptr);
258 
259 	switch (tag) {
260 	case TW_Denormal:
261 		if (denormal_operand() < 0)
262 			return;
263 		goto denormal_arg;
264 	case TW_Infinity:
265 		if (signnegative(st0_ptr)) {
266 			/* -infinity gives -1 (p16-10) */
267 			FPU_copy_to_reg0(&CONST_1, TAG_Valid);
268 			setnegative(st0_ptr);
269 		}
270 		return;
271 	default:
272 		single_arg_error(st0_ptr, tag);
273 	}
274 }
275 
276 static void fptan(FPU_REG *st0_ptr, u_char st0_tag)
277 {
278 	FPU_REG *st_new_ptr;
279 	int q;
280 	u_char arg_sign = getsign(st0_ptr);
281 
282 	/* Stack underflow has higher priority */
283 	if (st0_tag == TAG_Empty) {
284 		FPU_stack_underflow();	/* Puts a QNaN in st(0) */
285 		if (control_word & CW_Invalid) {
286 			st_new_ptr = &st(-1);
287 			push();
288 			FPU_stack_underflow();	/* Puts a QNaN in the new st(0) */
289 		}
290 		return;
291 	}
292 
293 	if (STACK_OVERFLOW) {
294 		FPU_stack_overflow();
295 		return;
296 	}
297 
298 	if (st0_tag == TAG_Valid) {
299 		if (exponent(st0_ptr) > -40) {
300 			if ((q = trig_arg(st0_ptr, 0)) == -1) {
301 				/* Operand is out of range */
302 				return;
303 			}
304 
305 			poly_tan(st0_ptr);
306 			setsign(st0_ptr, (q & 1) ^ (arg_sign != 0));
307 			set_precision_flag_up();	/* We do not really know if up or down */
308 		} else {
309 			/* For a small arg, the result == the argument */
310 			/* Underflow may happen */
311 
312 		      denormal_arg:
313 
314 			FPU_to_exp16(st0_ptr, st0_ptr);
315 
316 			st0_tag =
317 			    FPU_round(st0_ptr, 1, 0, FULL_PRECISION, arg_sign);
318 			FPU_settag0(st0_tag);
319 		}
320 		push();
321 		FPU_copy_to_reg0(&CONST_1, TAG_Valid);
322 		return;
323 	}
324 
325 	if (st0_tag == TAG_Zero) {
326 		push();
327 		FPU_copy_to_reg0(&CONST_1, TAG_Valid);
328 		setcc(0);
329 		return;
330 	}
331 
332 	if (st0_tag == TAG_Special)
333 		st0_tag = FPU_Special(st0_ptr);
334 
335 	if (st0_tag == TW_Denormal) {
336 		if (denormal_operand() < 0)
337 			return;
338 
339 		goto denormal_arg;
340 	}
341 
342 	if (st0_tag == TW_Infinity) {
343 		/* The 80486 treats infinity as an invalid operand */
344 		if (arith_invalid(0) >= 0) {
345 			st_new_ptr = &st(-1);
346 			push();
347 			arith_invalid(0);
348 		}
349 		return;
350 	}
351 
352 	single_arg_2_error(st0_ptr, st0_tag);
353 }
354 
355 static void fxtract(FPU_REG *st0_ptr, u_char st0_tag)
356 {
357 	FPU_REG *st_new_ptr;
358 	u_char sign;
359 	register FPU_REG *st1_ptr = st0_ptr;	/* anticipate */
360 
361 	if (STACK_OVERFLOW) {
362 		FPU_stack_overflow();
363 		return;
364 	}
365 
366 	clear_C1();
367 
368 	if (st0_tag == TAG_Valid) {
369 		long e;
370 
371 		push();
372 		sign = getsign(st1_ptr);
373 		reg_copy(st1_ptr, st_new_ptr);
374 		setexponent16(st_new_ptr, exponent(st_new_ptr));
375 
376 	      denormal_arg:
377 
378 		e = exponent16(st_new_ptr);
379 		convert_l2reg(&e, 1);
380 		setexponentpos(st_new_ptr, 0);
381 		setsign(st_new_ptr, sign);
382 		FPU_settag0(TAG_Valid);	/* Needed if arg was a denormal */
383 		return;
384 	} else if (st0_tag == TAG_Zero) {
385 		sign = getsign(st0_ptr);
386 
387 		if (FPU_divide_by_zero(0, SIGN_NEG) < 0)
388 			return;
389 
390 		push();
391 		FPU_copy_to_reg0(&CONST_Z, TAG_Zero);
392 		setsign(st_new_ptr, sign);
393 		return;
394 	}
395 
396 	if (st0_tag == TAG_Special)
397 		st0_tag = FPU_Special(st0_ptr);
398 
399 	if (st0_tag == TW_Denormal) {
400 		if (denormal_operand() < 0)
401 			return;
402 
403 		push();
404 		sign = getsign(st1_ptr);
405 		FPU_to_exp16(st1_ptr, st_new_ptr);
406 		goto denormal_arg;
407 	} else if (st0_tag == TW_Infinity) {
408 		sign = getsign(st0_ptr);
409 		setpositive(st0_ptr);
410 		push();
411 		FPU_copy_to_reg0(&CONST_INF, TAG_Special);
412 		setsign(st_new_ptr, sign);
413 		return;
414 	} else if (st0_tag == TW_NaN) {
415 		if (real_1op_NaN(st0_ptr) < 0)
416 			return;
417 
418 		push();
419 		FPU_copy_to_reg0(st0_ptr, TAG_Special);
420 		return;
421 	} else if (st0_tag == TAG_Empty) {
422 		/* Is this the correct behaviour? */
423 		if (control_word & EX_Invalid) {
424 			FPU_stack_underflow();
425 			push();
426 			FPU_stack_underflow();
427 		} else
428 			EXCEPTION(EX_StackUnder);
429 	}
430 #ifdef PARANOID
431 	else
432 		EXCEPTION(EX_INTERNAL | 0x119);
433 #endif /* PARANOID */
434 }
435 
436 static void fdecstp(void)
437 {
438 	clear_C1();
439 	top--;
440 }
441 
442 static void fincstp(void)
443 {
444 	clear_C1();
445 	top++;
446 }
447 
448 static void fsqrt_(FPU_REG *st0_ptr, u_char st0_tag)
449 {
450 	int expon;
451 
452 	clear_C1();
453 
454 	if (st0_tag == TAG_Valid) {
455 		u_char tag;
456 
457 		if (signnegative(st0_ptr)) {
458 			arith_invalid(0);	/* sqrt(negative) is invalid */
459 			return;
460 		}
461 
462 		/* make st(0) in  [1.0 .. 4.0) */
463 		expon = exponent(st0_ptr);
464 
465 	      denormal_arg:
466 
467 		setexponent16(st0_ptr, (expon & 1));
468 
469 		/* Do the computation, the sign of the result will be positive. */
470 		tag = wm_sqrt(st0_ptr, 0, 0, control_word, SIGN_POS);
471 		addexponent(st0_ptr, expon >> 1);
472 		FPU_settag0(tag);
473 		return;
474 	}
475 
476 	if (st0_tag == TAG_Zero)
477 		return;
478 
479 	if (st0_tag == TAG_Special)
480 		st0_tag = FPU_Special(st0_ptr);
481 
482 	if (st0_tag == TW_Infinity) {
483 		if (signnegative(st0_ptr))
484 			arith_invalid(0);	/* sqrt(-Infinity) is invalid */
485 		return;
486 	} else if (st0_tag == TW_Denormal) {
487 		if (signnegative(st0_ptr)) {
488 			arith_invalid(0);	/* sqrt(negative) is invalid */
489 			return;
490 		}
491 
492 		if (denormal_operand() < 0)
493 			return;
494 
495 		FPU_to_exp16(st0_ptr, st0_ptr);
496 
497 		expon = exponent16(st0_ptr);
498 
499 		goto denormal_arg;
500 	}
501 
502 	single_arg_error(st0_ptr, st0_tag);
503 
504 }
505 
506 static void frndint_(FPU_REG *st0_ptr, u_char st0_tag)
507 {
508 	int flags, tag;
509 
510 	if (st0_tag == TAG_Valid) {
511 		u_char sign;
512 
513 	      denormal_arg:
514 
515 		sign = getsign(st0_ptr);
516 
517 		if (exponent(st0_ptr) > 63)
518 			return;
519 
520 		if (st0_tag == TW_Denormal) {
521 			if (denormal_operand() < 0)
522 				return;
523 		}
524 
525 		/* Fortunately, this can't overflow to 2^64 */
526 		if ((flags = FPU_round_to_int(st0_ptr, st0_tag)))
527 			set_precision_flag(flags);
528 
529 		setexponent16(st0_ptr, 63);
530 		tag = FPU_normalize(st0_ptr);
531 		setsign(st0_ptr, sign);
532 		FPU_settag0(tag);
533 		return;
534 	}
535 
536 	if (st0_tag == TAG_Zero)
537 		return;
538 
539 	if (st0_tag == TAG_Special)
540 		st0_tag = FPU_Special(st0_ptr);
541 
542 	if (st0_tag == TW_Denormal)
543 		goto denormal_arg;
544 	else if (st0_tag == TW_Infinity)
545 		return;
546 	else
547 		single_arg_error(st0_ptr, st0_tag);
548 }
549 
550 static int fsin(FPU_REG *st0_ptr, u_char tag)
551 {
552 	u_char arg_sign = getsign(st0_ptr);
553 
554 	if (tag == TAG_Valid) {
555 		int q;
556 
557 		if (exponent(st0_ptr) > -40) {
558 			if ((q = trig_arg(st0_ptr, 0)) == -1) {
559 				/* Operand is out of range */
560 				return 1;
561 			}
562 
563 			poly_sine(st0_ptr);
564 
565 			if (q & 2)
566 				changesign(st0_ptr);
567 
568 			setsign(st0_ptr, getsign(st0_ptr) ^ arg_sign);
569 
570 			/* We do not really know if up or down */
571 			set_precision_flag_up();
572 			return 0;
573 		} else {
574 			/* For a small arg, the result == the argument */
575 			set_precision_flag_up();	/* Must be up. */
576 			return 0;
577 		}
578 	}
579 
580 	if (tag == TAG_Zero) {
581 		setcc(0);
582 		return 0;
583 	}
584 
585 	if (tag == TAG_Special)
586 		tag = FPU_Special(st0_ptr);
587 
588 	if (tag == TW_Denormal) {
589 		if (denormal_operand() < 0)
590 			return 1;
591 
592 		/* For a small arg, the result == the argument */
593 		/* Underflow may happen */
594 		FPU_to_exp16(st0_ptr, st0_ptr);
595 
596 		tag = FPU_round(st0_ptr, 1, 0, FULL_PRECISION, arg_sign);
597 
598 		FPU_settag0(tag);
599 
600 		return 0;
601 	} else if (tag == TW_Infinity) {
602 		/* The 80486 treats infinity as an invalid operand */
603 		arith_invalid(0);
604 		return 1;
605 	} else {
606 		single_arg_error(st0_ptr, tag);
607 		return 1;
608 	}
609 }
610 
611 static int f_cos(FPU_REG *st0_ptr, u_char tag)
612 {
613 	u_char st0_sign;
614 
615 	st0_sign = getsign(st0_ptr);
616 
617 	if (tag == TAG_Valid) {
618 		int q;
619 
620 		if (exponent(st0_ptr) > -40) {
621 			if ((exponent(st0_ptr) < 0)
622 			    || ((exponent(st0_ptr) == 0)
623 				&& (significand(st0_ptr) <=
624 				    0xc90fdaa22168c234LL))) {
625 				poly_cos(st0_ptr);
626 
627 				/* We do not really know if up or down */
628 				set_precision_flag_down();
629 
630 				return 0;
631 			} else if ((q = trig_arg(st0_ptr, FCOS)) != -1) {
632 				poly_sine(st0_ptr);
633 
634 				if ((q + 1) & 2)
635 					changesign(st0_ptr);
636 
637 				/* We do not really know if up or down */
638 				set_precision_flag_down();
639 
640 				return 0;
641 			} else {
642 				/* Operand is out of range */
643 				return 1;
644 			}
645 		} else {
646 		      denormal_arg:
647 
648 			setcc(0);
649 			FPU_copy_to_reg0(&CONST_1, TAG_Valid);
650 #ifdef PECULIAR_486
651 			set_precision_flag_down();	/* 80486 appears to do this. */
652 #else
653 			set_precision_flag_up();	/* Must be up. */
654 #endif /* PECULIAR_486 */
655 			return 0;
656 		}
657 	} else if (tag == TAG_Zero) {
658 		FPU_copy_to_reg0(&CONST_1, TAG_Valid);
659 		setcc(0);
660 		return 0;
661 	}
662 
663 	if (tag == TAG_Special)
664 		tag = FPU_Special(st0_ptr);
665 
666 	if (tag == TW_Denormal) {
667 		if (denormal_operand() < 0)
668 			return 1;
669 
670 		goto denormal_arg;
671 	} else if (tag == TW_Infinity) {
672 		/* The 80486 treats infinity as an invalid operand */
673 		arith_invalid(0);
674 		return 1;
675 	} else {
676 		single_arg_error(st0_ptr, tag);	/* requires st0_ptr == &st(0) */
677 		return 1;
678 	}
679 }
680 
681 static void fcos(FPU_REG *st0_ptr, u_char st0_tag)
682 {
683 	f_cos(st0_ptr, st0_tag);
684 }
685 
686 static void fsincos(FPU_REG *st0_ptr, u_char st0_tag)
687 {
688 	FPU_REG *st_new_ptr;
689 	FPU_REG arg;
690 	u_char tag;
691 
692 	/* Stack underflow has higher priority */
693 	if (st0_tag == TAG_Empty) {
694 		FPU_stack_underflow();	/* Puts a QNaN in st(0) */
695 		if (control_word & CW_Invalid) {
696 			st_new_ptr = &st(-1);
697 			push();
698 			FPU_stack_underflow();	/* Puts a QNaN in the new st(0) */
699 		}
700 		return;
701 	}
702 
703 	if (STACK_OVERFLOW) {
704 		FPU_stack_overflow();
705 		return;
706 	}
707 
708 	if (st0_tag == TAG_Special)
709 		tag = FPU_Special(st0_ptr);
710 	else
711 		tag = st0_tag;
712 
713 	if (tag == TW_NaN) {
714 		single_arg_2_error(st0_ptr, TW_NaN);
715 		return;
716 	} else if (tag == TW_Infinity) {
717 		/* The 80486 treats infinity as an invalid operand */
718 		if (arith_invalid(0) >= 0) {
719 			/* Masked response */
720 			push();
721 			arith_invalid(0);
722 		}
723 		return;
724 	}
725 
726 	reg_copy(st0_ptr, &arg);
727 	if (!fsin(st0_ptr, st0_tag)) {
728 		push();
729 		FPU_copy_to_reg0(&arg, st0_tag);
730 		f_cos(&st(0), st0_tag);
731 	} else {
732 		/* An error, so restore st(0) */
733 		FPU_copy_to_reg0(&arg, st0_tag);
734 	}
735 }
736 
737 /*---------------------------------------------------------------------------*/
738 /* The following all require two arguments: st(0) and st(1) */
739 
740 /* A lean, mean kernel for the fprem instructions. This relies upon
741    the division and rounding to an integer in do_fprem giving an
742    exact result. Because of this, rem_kernel() needs to deal only with
743    the least significant 64 bits, the more significant bits of the
744    result must be zero.
745  */
746 static void rem_kernel(unsigned long long st0, unsigned long long *y,
747 		       unsigned long long st1, unsigned long long q, int n)
748 {
749 	int dummy;
750 	unsigned long long x;
751 
752 	x = st0 << n;
753 
754 	/* Do the required multiplication and subtraction in the one operation */
755 
756 	/* lsw x -= lsw st1 * lsw q */
757 	asm volatile ("mull %4; subl %%eax,%0; sbbl %%edx,%1":"=m"
758 		      (((unsigned *)&x)[0]), "=m"(((unsigned *)&x)[1]),
759 		      "=a"(dummy)
760 		      :"2"(((unsigned *)&st1)[0]), "m"(((unsigned *)&q)[0])
761 		      :"%dx");
762 	/* msw x -= msw st1 * lsw q */
763 	asm volatile ("mull %3; subl %%eax,%0":"=m" (((unsigned *)&x)[1]),
764 		      "=a"(dummy)
765 		      :"1"(((unsigned *)&st1)[1]), "m"(((unsigned *)&q)[0])
766 		      :"%dx");
767 	/* msw x -= lsw st1 * msw q */
768 	asm volatile ("mull %3; subl %%eax,%0":"=m" (((unsigned *)&x)[1]),
769 		      "=a"(dummy)
770 		      :"1"(((unsigned *)&st1)[0]), "m"(((unsigned *)&q)[1])
771 		      :"%dx");
772 
773 	*y = x;
774 }
775 
776 /* Remainder of st(0) / st(1) */
777 /* This routine produces exact results, i.e. there is never any
778    rounding or truncation, etc of the result. */
779 static void do_fprem(FPU_REG *st0_ptr, u_char st0_tag, int round)
780 {
781 	FPU_REG *st1_ptr = &st(1);
782 	u_char st1_tag = FPU_gettagi(1);
783 
784 	if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
785 		FPU_REG tmp, st0, st1;
786 		u_char st0_sign, st1_sign;
787 		u_char tmptag;
788 		int tag;
789 		int old_cw;
790 		int expdif;
791 		long long q;
792 		unsigned short saved_status;
793 		int cc;
794 
795 	      fprem_valid:
796 		/* Convert registers for internal use. */
797 		st0_sign = FPU_to_exp16(st0_ptr, &st0);
798 		st1_sign = FPU_to_exp16(st1_ptr, &st1);
799 		expdif = exponent16(&st0) - exponent16(&st1);
800 
801 		old_cw = control_word;
802 		cc = 0;
803 
804 		/* We want the status following the denorm tests, but don't want
805 		   the status changed by the arithmetic operations. */
806 		saved_status = partial_status;
807 		control_word &= ~CW_RC;
808 		control_word |= RC_CHOP;
809 
810 		if (expdif < 64) {
811 			/* This should be the most common case */
812 
813 			if (expdif > -2) {
814 				u_char sign = st0_sign ^ st1_sign;
815 				tag = FPU_u_div(&st0, &st1, &tmp,
816 						PR_64_BITS | RC_CHOP | 0x3f,
817 						sign);
818 				setsign(&tmp, sign);
819 
820 				if (exponent(&tmp) >= 0) {
821 					FPU_round_to_int(&tmp, tag);	/* Fortunately, this can't
822 									   overflow to 2^64 */
823 					q = significand(&tmp);
824 
825 					rem_kernel(significand(&st0),
826 						   &significand(&tmp),
827 						   significand(&st1),
828 						   q, expdif);
829 
830 					setexponent16(&tmp, exponent16(&st1));
831 				} else {
832 					reg_copy(&st0, &tmp);
833 					q = 0;
834 				}
835 
836 				if ((round == RC_RND)
837 				    && (tmp.sigh & 0xc0000000)) {
838 					/* We may need to subtract st(1) once more,
839 					   to get a result <= 1/2 of st(1). */
840 					unsigned long long x;
841 					expdif =
842 					    exponent16(&st1) - exponent16(&tmp);
843 					if (expdif <= 1) {
844 						if (expdif == 0)
845 							x = significand(&st1) -
846 							    significand(&tmp);
847 						else	/* expdif is 1 */
848 							x = (significand(&st1)
849 							     << 1) -
850 							    significand(&tmp);
851 						if ((x < significand(&tmp)) ||
852 						    /* or equi-distant (from 0 & st(1)) and q is odd */
853 						    ((x == significand(&tmp))
854 						     && (q & 1))) {
855 							st0_sign = !st0_sign;
856 							significand(&tmp) = x;
857 							q++;
858 						}
859 					}
860 				}
861 
862 				if (q & 4)
863 					cc |= SW_C0;
864 				if (q & 2)
865 					cc |= SW_C3;
866 				if (q & 1)
867 					cc |= SW_C1;
868 			} else {
869 				control_word = old_cw;
870 				setcc(0);
871 				return;
872 			}
873 		} else {
874 			/* There is a large exponent difference ( >= 64 ) */
875 			/* To make much sense, the code in this section should
876 			   be done at high precision. */
877 			int exp_1, N;
878 			u_char sign;
879 
880 			/* prevent overflow here */
881 			/* N is 'a number between 32 and 63' (p26-113) */
882 			reg_copy(&st0, &tmp);
883 			tmptag = st0_tag;
884 			N = (expdif & 0x0000001f) + 32;	/* This choice gives results
885 							   identical to an AMD 486 */
886 			setexponent16(&tmp, N);
887 			exp_1 = exponent16(&st1);
888 			setexponent16(&st1, 0);
889 			expdif -= N;
890 
891 			sign = getsign(&tmp) ^ st1_sign;
892 			tag =
893 			    FPU_u_div(&tmp, &st1, &tmp,
894 				      PR_64_BITS | RC_CHOP | 0x3f, sign);
895 			setsign(&tmp, sign);
896 
897 			FPU_round_to_int(&tmp, tag);	/* Fortunately, this can't
898 							   overflow to 2^64 */
899 
900 			rem_kernel(significand(&st0),
901 				   &significand(&tmp),
902 				   significand(&st1),
903 				   significand(&tmp), exponent(&tmp)
904 			    );
905 			setexponent16(&tmp, exp_1 + expdif);
906 
907 			/* It is possible for the operation to be complete here.
908 			   What does the IEEE standard say? The Intel 80486 manual
909 			   implies that the operation will never be completed at this
910 			   point, and the behaviour of a real 80486 confirms this.
911 			 */
912 			if (!(tmp.sigh | tmp.sigl)) {
913 				/* The result is zero */
914 				control_word = old_cw;
915 				partial_status = saved_status;
916 				FPU_copy_to_reg0(&CONST_Z, TAG_Zero);
917 				setsign(&st0, st0_sign);
918 #ifdef PECULIAR_486
919 				setcc(SW_C2);
920 #else
921 				setcc(0);
922 #endif /* PECULIAR_486 */
923 				return;
924 			}
925 			cc = SW_C2;
926 		}
927 
928 		control_word = old_cw;
929 		partial_status = saved_status;
930 		tag = FPU_normalize_nuo(&tmp);
931 		reg_copy(&tmp, st0_ptr);
932 
933 		/* The only condition to be looked for is underflow,
934 		   and it can occur here only if underflow is unmasked. */
935 		if ((exponent16(&tmp) <= EXP_UNDER) && (tag != TAG_Zero)
936 		    && !(control_word & CW_Underflow)) {
937 			setcc(cc);
938 			tag = arith_underflow(st0_ptr);
939 			setsign(st0_ptr, st0_sign);
940 			FPU_settag0(tag);
941 			return;
942 		} else if ((exponent16(&tmp) > EXP_UNDER) || (tag == TAG_Zero)) {
943 			stdexp(st0_ptr);
944 			setsign(st0_ptr, st0_sign);
945 		} else {
946 			tag =
947 			    FPU_round(st0_ptr, 0, 0, FULL_PRECISION, st0_sign);
948 		}
949 		FPU_settag0(tag);
950 		setcc(cc);
951 
952 		return;
953 	}
954 
955 	if (st0_tag == TAG_Special)
956 		st0_tag = FPU_Special(st0_ptr);
957 	if (st1_tag == TAG_Special)
958 		st1_tag = FPU_Special(st1_ptr);
959 
960 	if (((st0_tag == TAG_Valid) && (st1_tag == TW_Denormal))
961 	    || ((st0_tag == TW_Denormal) && (st1_tag == TAG_Valid))
962 	    || ((st0_tag == TW_Denormal) && (st1_tag == TW_Denormal))) {
963 		if (denormal_operand() < 0)
964 			return;
965 		goto fprem_valid;
966 	} else if ((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty)) {
967 		FPU_stack_underflow();
968 		return;
969 	} else if (st0_tag == TAG_Zero) {
970 		if (st1_tag == TAG_Valid) {
971 			setcc(0);
972 			return;
973 		} else if (st1_tag == TW_Denormal) {
974 			if (denormal_operand() < 0)
975 				return;
976 			setcc(0);
977 			return;
978 		} else if (st1_tag == TAG_Zero) {
979 			arith_invalid(0);
980 			return;
981 		} /* fprem(?,0) always invalid */
982 		else if (st1_tag == TW_Infinity) {
983 			setcc(0);
984 			return;
985 		}
986 	} else if ((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal)) {
987 		if (st1_tag == TAG_Zero) {
988 			arith_invalid(0);	/* fprem(Valid,Zero) is invalid */
989 			return;
990 		} else if (st1_tag != TW_NaN) {
991 			if (((st0_tag == TW_Denormal)
992 			     || (st1_tag == TW_Denormal))
993 			    && (denormal_operand() < 0))
994 				return;
995 
996 			if (st1_tag == TW_Infinity) {
997 				/* fprem(Valid,Infinity) is o.k. */
998 				setcc(0);
999 				return;
1000 			}
1001 		}
1002 	} else if (st0_tag == TW_Infinity) {
1003 		if (st1_tag != TW_NaN) {
1004 			arith_invalid(0);	/* fprem(Infinity,?) is invalid */
1005 			return;
1006 		}
1007 	}
1008 
1009 	/* One of the registers must contain a NaN if we got here. */
1010 
1011 #ifdef PARANOID
1012 	if ((st0_tag != TW_NaN) && (st1_tag != TW_NaN))
1013 		EXCEPTION(EX_INTERNAL | 0x118);
1014 #endif /* PARANOID */
1015 
1016 	real_2op_NaN(st1_ptr, st1_tag, 0, st1_ptr);
1017 
1018 }
1019 
1020 /* ST(1) <- ST(1) * log ST;  pop ST */
1021 static void fyl2x(FPU_REG *st0_ptr, u_char st0_tag)
1022 {
1023 	FPU_REG *st1_ptr = &st(1), exponent;
1024 	u_char st1_tag = FPU_gettagi(1);
1025 	u_char sign;
1026 	int e, tag;
1027 
1028 	clear_C1();
1029 
1030 	if ((st0_tag == TAG_Valid) && (st1_tag == TAG_Valid)) {
1031 	      both_valid:
1032 		/* Both regs are Valid or Denormal */
1033 		if (signpositive(st0_ptr)) {
1034 			if (st0_tag == TW_Denormal)
1035 				FPU_to_exp16(st0_ptr, st0_ptr);
1036 			else
1037 				/* Convert st(0) for internal use. */
1038 				setexponent16(st0_ptr, exponent(st0_ptr));
1039 
1040 			if ((st0_ptr->sigh == 0x80000000)
1041 			    && (st0_ptr->sigl == 0)) {
1042 				/* Special case. The result can be precise. */
1043 				u_char esign;
1044 				e = exponent16(st0_ptr);
1045 				if (e >= 0) {
1046 					exponent.sigh = e;
1047 					esign = SIGN_POS;
1048 				} else {
1049 					exponent.sigh = -e;
1050 					esign = SIGN_NEG;
1051 				}
1052 				exponent.sigl = 0;
1053 				setexponent16(&exponent, 31);
1054 				tag = FPU_normalize_nuo(&exponent);
1055 				stdexp(&exponent);
1056 				setsign(&exponent, esign);
1057 				tag =
1058 				    FPU_mul(&exponent, tag, 1, FULL_PRECISION);
1059 				if (tag >= 0)
1060 					FPU_settagi(1, tag);
1061 			} else {
1062 				/* The usual case */
1063 				sign = getsign(st1_ptr);
1064 				if (st1_tag == TW_Denormal)
1065 					FPU_to_exp16(st1_ptr, st1_ptr);
1066 				else
1067 					/* Convert st(1) for internal use. */
1068 					setexponent16(st1_ptr,
1069 						      exponent(st1_ptr));
1070 				poly_l2(st0_ptr, st1_ptr, sign);
1071 			}
1072 		} else {
1073 			/* negative */
1074 			if (arith_invalid(1) < 0)
1075 				return;
1076 		}
1077 
1078 		FPU_pop();
1079 
1080 		return;
1081 	}
1082 
1083 	if (st0_tag == TAG_Special)
1084 		st0_tag = FPU_Special(st0_ptr);
1085 	if (st1_tag == TAG_Special)
1086 		st1_tag = FPU_Special(st1_ptr);
1087 
1088 	if ((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty)) {
1089 		FPU_stack_underflow_pop(1);
1090 		return;
1091 	} else if ((st0_tag <= TW_Denormal) && (st1_tag <= TW_Denormal)) {
1092 		if (st0_tag == TAG_Zero) {
1093 			if (st1_tag == TAG_Zero) {
1094 				/* Both args zero is invalid */
1095 				if (arith_invalid(1) < 0)
1096 					return;
1097 			} else {
1098 				u_char sign;
1099 				sign = getsign(st1_ptr) ^ SIGN_NEG;
1100 				if (FPU_divide_by_zero(1, sign) < 0)
1101 					return;
1102 
1103 				setsign(st1_ptr, sign);
1104 			}
1105 		} else if (st1_tag == TAG_Zero) {
1106 			/* st(1) contains zero, st(0) valid <> 0 */
1107 			/* Zero is the valid answer */
1108 			sign = getsign(st1_ptr);
1109 
1110 			if (signnegative(st0_ptr)) {
1111 				/* log(negative) */
1112 				if (arith_invalid(1) < 0)
1113 					return;
1114 			} else if ((st0_tag == TW_Denormal)
1115 				   && (denormal_operand() < 0))
1116 				return;
1117 			else {
1118 				if (exponent(st0_ptr) < 0)
1119 					sign ^= SIGN_NEG;
1120 
1121 				FPU_copy_to_reg1(&CONST_Z, TAG_Zero);
1122 				setsign(st1_ptr, sign);
1123 			}
1124 		} else {
1125 			/* One or both operands are denormals. */
1126 			if (denormal_operand() < 0)
1127 				return;
1128 			goto both_valid;
1129 		}
1130 	} else if ((st0_tag == TW_NaN) || (st1_tag == TW_NaN)) {
1131 		if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1132 			return;
1133 	}
1134 	/* One or both arg must be an infinity */
1135 	else if (st0_tag == TW_Infinity) {
1136 		if ((signnegative(st0_ptr)) || (st1_tag == TAG_Zero)) {
1137 			/* log(-infinity) or 0*log(infinity) */
1138 			if (arith_invalid(1) < 0)
1139 				return;
1140 		} else {
1141 			u_char sign = getsign(st1_ptr);
1142 
1143 			if ((st1_tag == TW_Denormal)
1144 			    && (denormal_operand() < 0))
1145 				return;
1146 
1147 			FPU_copy_to_reg1(&CONST_INF, TAG_Special);
1148 			setsign(st1_ptr, sign);
1149 		}
1150 	}
1151 	/* st(1) must be infinity here */
1152 	else if (((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal))
1153 		 && (signpositive(st0_ptr))) {
1154 		if (exponent(st0_ptr) >= 0) {
1155 			if ((exponent(st0_ptr) == 0) &&
1156 			    (st0_ptr->sigh == 0x80000000) &&
1157 			    (st0_ptr->sigl == 0)) {
1158 				/* st(0) holds 1.0 */
1159 				/* infinity*log(1) */
1160 				if (arith_invalid(1) < 0)
1161 					return;
1162 			}
1163 			/* else st(0) is positive and > 1.0 */
1164 		} else {
1165 			/* st(0) is positive and < 1.0 */
1166 
1167 			if ((st0_tag == TW_Denormal)
1168 			    && (denormal_operand() < 0))
1169 				return;
1170 
1171 			changesign(st1_ptr);
1172 		}
1173 	} else {
1174 		/* st(0) must be zero or negative */
1175 		if (st0_tag == TAG_Zero) {
1176 			/* This should be invalid, but a real 80486 is happy with it. */
1177 
1178 #ifndef PECULIAR_486
1179 			sign = getsign(st1_ptr);
1180 			if (FPU_divide_by_zero(1, sign) < 0)
1181 				return;
1182 #endif /* PECULIAR_486 */
1183 
1184 			changesign(st1_ptr);
1185 		} else if (arith_invalid(1) < 0)	/* log(negative) */
1186 			return;
1187 	}
1188 
1189 	FPU_pop();
1190 }
1191 
1192 static void fpatan(FPU_REG *st0_ptr, u_char st0_tag)
1193 {
1194 	FPU_REG *st1_ptr = &st(1);
1195 	u_char st1_tag = FPU_gettagi(1);
1196 	int tag;
1197 
1198 	clear_C1();
1199 	if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
1200 	      valid_atan:
1201 
1202 		poly_atan(st0_ptr, st0_tag, st1_ptr, st1_tag);
1203 
1204 		FPU_pop();
1205 
1206 		return;
1207 	}
1208 
1209 	if (st0_tag == TAG_Special)
1210 		st0_tag = FPU_Special(st0_ptr);
1211 	if (st1_tag == TAG_Special)
1212 		st1_tag = FPU_Special(st1_ptr);
1213 
1214 	if (((st0_tag == TAG_Valid) && (st1_tag == TW_Denormal))
1215 	    || ((st0_tag == TW_Denormal) && (st1_tag == TAG_Valid))
1216 	    || ((st0_tag == TW_Denormal) && (st1_tag == TW_Denormal))) {
1217 		if (denormal_operand() < 0)
1218 			return;
1219 
1220 		goto valid_atan;
1221 	} else if ((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty)) {
1222 		FPU_stack_underflow_pop(1);
1223 		return;
1224 	} else if ((st0_tag == TW_NaN) || (st1_tag == TW_NaN)) {
1225 		if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) >= 0)
1226 			FPU_pop();
1227 		return;
1228 	} else if ((st0_tag == TW_Infinity) || (st1_tag == TW_Infinity)) {
1229 		u_char sign = getsign(st1_ptr);
1230 		if (st0_tag == TW_Infinity) {
1231 			if (st1_tag == TW_Infinity) {
1232 				if (signpositive(st0_ptr)) {
1233 					FPU_copy_to_reg1(&CONST_PI4, TAG_Valid);
1234 				} else {
1235 					setpositive(st1_ptr);
1236 					tag =
1237 					    FPU_u_add(&CONST_PI4, &CONST_PI2,
1238 						      st1_ptr, FULL_PRECISION,
1239 						      SIGN_POS,
1240 						      exponent(&CONST_PI4),
1241 						      exponent(&CONST_PI2));
1242 					if (tag >= 0)
1243 						FPU_settagi(1, tag);
1244 				}
1245 			} else {
1246 				if ((st1_tag == TW_Denormal)
1247 				    && (denormal_operand() < 0))
1248 					return;
1249 
1250 				if (signpositive(st0_ptr)) {
1251 					FPU_copy_to_reg1(&CONST_Z, TAG_Zero);
1252 					setsign(st1_ptr, sign);	/* An 80486 preserves the sign */
1253 					FPU_pop();
1254 					return;
1255 				} else {
1256 					FPU_copy_to_reg1(&CONST_PI, TAG_Valid);
1257 				}
1258 			}
1259 		} else {
1260 			/* st(1) is infinity, st(0) not infinity */
1261 			if ((st0_tag == TW_Denormal)
1262 			    && (denormal_operand() < 0))
1263 				return;
1264 
1265 			FPU_copy_to_reg1(&CONST_PI2, TAG_Valid);
1266 		}
1267 		setsign(st1_ptr, sign);
1268 	} else if (st1_tag == TAG_Zero) {
1269 		/* st(0) must be valid or zero */
1270 		u_char sign = getsign(st1_ptr);
1271 
1272 		if ((st0_tag == TW_Denormal) && (denormal_operand() < 0))
1273 			return;
1274 
1275 		if (signpositive(st0_ptr)) {
1276 			/* An 80486 preserves the sign */
1277 			FPU_pop();
1278 			return;
1279 		}
1280 
1281 		FPU_copy_to_reg1(&CONST_PI, TAG_Valid);
1282 		setsign(st1_ptr, sign);
1283 	} else if (st0_tag == TAG_Zero) {
1284 		/* st(1) must be TAG_Valid here */
1285 		u_char sign = getsign(st1_ptr);
1286 
1287 		if ((st1_tag == TW_Denormal) && (denormal_operand() < 0))
1288 			return;
1289 
1290 		FPU_copy_to_reg1(&CONST_PI2, TAG_Valid);
1291 		setsign(st1_ptr, sign);
1292 	}
1293 #ifdef PARANOID
1294 	else
1295 		EXCEPTION(EX_INTERNAL | 0x125);
1296 #endif /* PARANOID */
1297 
1298 	FPU_pop();
1299 	set_precision_flag_up();	/* We do not really know if up or down */
1300 }
1301 
1302 static void fprem(FPU_REG *st0_ptr, u_char st0_tag)
1303 {
1304 	do_fprem(st0_ptr, st0_tag, RC_CHOP);
1305 }
1306 
1307 static void fprem1(FPU_REG *st0_ptr, u_char st0_tag)
1308 {
1309 	do_fprem(st0_ptr, st0_tag, RC_RND);
1310 }
1311 
1312 static void fyl2xp1(FPU_REG *st0_ptr, u_char st0_tag)
1313 {
1314 	u_char sign, sign1;
1315 	FPU_REG *st1_ptr = &st(1), a, b;
1316 	u_char st1_tag = FPU_gettagi(1);
1317 
1318 	clear_C1();
1319 	if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
1320 	      valid_yl2xp1:
1321 
1322 		sign = getsign(st0_ptr);
1323 		sign1 = getsign(st1_ptr);
1324 
1325 		FPU_to_exp16(st0_ptr, &a);
1326 		FPU_to_exp16(st1_ptr, &b);
1327 
1328 		if (poly_l2p1(sign, sign1, &a, &b, st1_ptr))
1329 			return;
1330 
1331 		FPU_pop();
1332 		return;
1333 	}
1334 
1335 	if (st0_tag == TAG_Special)
1336 		st0_tag = FPU_Special(st0_ptr);
1337 	if (st1_tag == TAG_Special)
1338 		st1_tag = FPU_Special(st1_ptr);
1339 
1340 	if (((st0_tag == TAG_Valid) && (st1_tag == TW_Denormal))
1341 	    || ((st0_tag == TW_Denormal) && (st1_tag == TAG_Valid))
1342 	    || ((st0_tag == TW_Denormal) && (st1_tag == TW_Denormal))) {
1343 		if (denormal_operand() < 0)
1344 			return;
1345 
1346 		goto valid_yl2xp1;
1347 	} else if ((st0_tag == TAG_Empty) | (st1_tag == TAG_Empty)) {
1348 		FPU_stack_underflow_pop(1);
1349 		return;
1350 	} else if (st0_tag == TAG_Zero) {
1351 		switch (st1_tag) {
1352 		case TW_Denormal:
1353 			if (denormal_operand() < 0)
1354 				return;
1355 			fallthrough;
1356 		case TAG_Zero:
1357 		case TAG_Valid:
1358 			setsign(st0_ptr, getsign(st0_ptr) ^ getsign(st1_ptr));
1359 			FPU_copy_to_reg1(st0_ptr, st0_tag);
1360 			break;
1361 
1362 		case TW_Infinity:
1363 			/* Infinity*log(1) */
1364 			if (arith_invalid(1) < 0)
1365 				return;
1366 			break;
1367 
1368 		case TW_NaN:
1369 			if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1370 				return;
1371 			break;
1372 
1373 		default:
1374 #ifdef PARANOID
1375 			EXCEPTION(EX_INTERNAL | 0x116);
1376 			return;
1377 #endif /* PARANOID */
1378 			break;
1379 		}
1380 	} else if ((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal)) {
1381 		switch (st1_tag) {
1382 		case TAG_Zero:
1383 			if (signnegative(st0_ptr)) {
1384 				if (exponent(st0_ptr) >= 0) {
1385 					/* st(0) holds <= -1.0 */
1386 #ifdef PECULIAR_486		/* Stupid 80486 doesn't worry about log(negative). */
1387 					changesign(st1_ptr);
1388 #else
1389 					if (arith_invalid(1) < 0)
1390 						return;
1391 #endif /* PECULIAR_486 */
1392 				} else if ((st0_tag == TW_Denormal)
1393 					   && (denormal_operand() < 0))
1394 					return;
1395 				else
1396 					changesign(st1_ptr);
1397 			} else if ((st0_tag == TW_Denormal)
1398 				   && (denormal_operand() < 0))
1399 				return;
1400 			break;
1401 
1402 		case TW_Infinity:
1403 			if (signnegative(st0_ptr)) {
1404 				if ((exponent(st0_ptr) >= 0) &&
1405 				    !((st0_ptr->sigh == 0x80000000) &&
1406 				      (st0_ptr->sigl == 0))) {
1407 					/* st(0) holds < -1.0 */
1408 #ifdef PECULIAR_486		/* Stupid 80486 doesn't worry about log(negative). */
1409 					changesign(st1_ptr);
1410 #else
1411 					if (arith_invalid(1) < 0)
1412 						return;
1413 #endif /* PECULIAR_486 */
1414 				} else if ((st0_tag == TW_Denormal)
1415 					   && (denormal_operand() < 0))
1416 					return;
1417 				else
1418 					changesign(st1_ptr);
1419 			} else if ((st0_tag == TW_Denormal)
1420 				   && (denormal_operand() < 0))
1421 				return;
1422 			break;
1423 
1424 		case TW_NaN:
1425 			if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1426 				return;
1427 		}
1428 
1429 	} else if (st0_tag == TW_NaN) {
1430 		if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1431 			return;
1432 	} else if (st0_tag == TW_Infinity) {
1433 		if (st1_tag == TW_NaN) {
1434 			if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1435 				return;
1436 		} else if (signnegative(st0_ptr)) {
1437 #ifndef PECULIAR_486
1438 			/* This should have higher priority than denormals, but... */
1439 			if (arith_invalid(1) < 0)	/* log(-infinity) */
1440 				return;
1441 #endif /* PECULIAR_486 */
1442 			if ((st1_tag == TW_Denormal)
1443 			    && (denormal_operand() < 0))
1444 				return;
1445 #ifdef PECULIAR_486
1446 			/* Denormal operands actually get higher priority */
1447 			if (arith_invalid(1) < 0)	/* log(-infinity) */
1448 				return;
1449 #endif /* PECULIAR_486 */
1450 		} else if (st1_tag == TAG_Zero) {
1451 			/* log(infinity) */
1452 			if (arith_invalid(1) < 0)
1453 				return;
1454 		}
1455 
1456 		/* st(1) must be valid here. */
1457 
1458 		else if ((st1_tag == TW_Denormal) && (denormal_operand() < 0))
1459 			return;
1460 
1461 		/* The Manual says that log(Infinity) is invalid, but a real
1462 		   80486 sensibly says that it is o.k. */
1463 		else {
1464 			u_char sign = getsign(st1_ptr);
1465 			FPU_copy_to_reg1(&CONST_INF, TAG_Special);
1466 			setsign(st1_ptr, sign);
1467 		}
1468 	}
1469 #ifdef PARANOID
1470 	else {
1471 		EXCEPTION(EX_INTERNAL | 0x117);
1472 		return;
1473 	}
1474 #endif /* PARANOID */
1475 
1476 	FPU_pop();
1477 	return;
1478 
1479 }
1480 
1481 static void fscale(FPU_REG *st0_ptr, u_char st0_tag)
1482 {
1483 	FPU_REG *st1_ptr = &st(1);
1484 	u_char st1_tag = FPU_gettagi(1);
1485 	int old_cw = control_word;
1486 	u_char sign = getsign(st0_ptr);
1487 
1488 	clear_C1();
1489 	if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
1490 		long scale;
1491 		FPU_REG tmp;
1492 
1493 		/* Convert register for internal use. */
1494 		setexponent16(st0_ptr, exponent(st0_ptr));
1495 
1496 	      valid_scale:
1497 
1498 		if (exponent(st1_ptr) > 30) {
1499 			/* 2^31 is far too large, would require 2^(2^30) or 2^(-2^30) */
1500 
1501 			if (signpositive(st1_ptr)) {
1502 				EXCEPTION(EX_Overflow);
1503 				FPU_copy_to_reg0(&CONST_INF, TAG_Special);
1504 			} else {
1505 				EXCEPTION(EX_Underflow);
1506 				FPU_copy_to_reg0(&CONST_Z, TAG_Zero);
1507 			}
1508 			setsign(st0_ptr, sign);
1509 			return;
1510 		}
1511 
1512 		control_word &= ~CW_RC;
1513 		control_word |= RC_CHOP;
1514 		reg_copy(st1_ptr, &tmp);
1515 		FPU_round_to_int(&tmp, st1_tag);	/* This can never overflow here */
1516 		control_word = old_cw;
1517 		scale = signnegative(st1_ptr) ? -tmp.sigl : tmp.sigl;
1518 		scale += exponent16(st0_ptr);
1519 
1520 		setexponent16(st0_ptr, scale);
1521 
1522 		/* Use FPU_round() to properly detect under/overflow etc */
1523 		FPU_round(st0_ptr, 0, 0, control_word, sign);
1524 
1525 		return;
1526 	}
1527 
1528 	if (st0_tag == TAG_Special)
1529 		st0_tag = FPU_Special(st0_ptr);
1530 	if (st1_tag == TAG_Special)
1531 		st1_tag = FPU_Special(st1_ptr);
1532 
1533 	if ((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal)) {
1534 		switch (st1_tag) {
1535 		case TAG_Valid:
1536 			/* st(0) must be a denormal */
1537 			if ((st0_tag == TW_Denormal)
1538 			    && (denormal_operand() < 0))
1539 				return;
1540 
1541 			FPU_to_exp16(st0_ptr, st0_ptr);	/* Will not be left on stack */
1542 			goto valid_scale;
1543 
1544 		case TAG_Zero:
1545 			if (st0_tag == TW_Denormal)
1546 				denormal_operand();
1547 			return;
1548 
1549 		case TW_Denormal:
1550 			denormal_operand();
1551 			return;
1552 
1553 		case TW_Infinity:
1554 			if ((st0_tag == TW_Denormal)
1555 			    && (denormal_operand() < 0))
1556 				return;
1557 
1558 			if (signpositive(st1_ptr))
1559 				FPU_copy_to_reg0(&CONST_INF, TAG_Special);
1560 			else
1561 				FPU_copy_to_reg0(&CONST_Z, TAG_Zero);
1562 			setsign(st0_ptr, sign);
1563 			return;
1564 
1565 		case TW_NaN:
1566 			real_2op_NaN(st1_ptr, st1_tag, 0, st0_ptr);
1567 			return;
1568 		}
1569 	} else if (st0_tag == TAG_Zero) {
1570 		switch (st1_tag) {
1571 		case TAG_Valid:
1572 		case TAG_Zero:
1573 			return;
1574 
1575 		case TW_Denormal:
1576 			denormal_operand();
1577 			return;
1578 
1579 		case TW_Infinity:
1580 			if (signpositive(st1_ptr))
1581 				arith_invalid(0);	/* Zero scaled by +Infinity */
1582 			return;
1583 
1584 		case TW_NaN:
1585 			real_2op_NaN(st1_ptr, st1_tag, 0, st0_ptr);
1586 			return;
1587 		}
1588 	} else if (st0_tag == TW_Infinity) {
1589 		switch (st1_tag) {
1590 		case TAG_Valid:
1591 		case TAG_Zero:
1592 			return;
1593 
1594 		case TW_Denormal:
1595 			denormal_operand();
1596 			return;
1597 
1598 		case TW_Infinity:
1599 			if (signnegative(st1_ptr))
1600 				arith_invalid(0);	/* Infinity scaled by -Infinity */
1601 			return;
1602 
1603 		case TW_NaN:
1604 			real_2op_NaN(st1_ptr, st1_tag, 0, st0_ptr);
1605 			return;
1606 		}
1607 	} else if (st0_tag == TW_NaN) {
1608 		if (st1_tag != TAG_Empty) {
1609 			real_2op_NaN(st1_ptr, st1_tag, 0, st0_ptr);
1610 			return;
1611 		}
1612 	}
1613 #ifdef PARANOID
1614 	if (!((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty))) {
1615 		EXCEPTION(EX_INTERNAL | 0x115);
1616 		return;
1617 	}
1618 #endif
1619 
1620 	/* At least one of st(0), st(1) must be empty */
1621 	FPU_stack_underflow();
1622 
1623 }
1624 
1625 /*---------------------------------------------------------------------------*/
1626 
1627 static FUNC_ST0 const trig_table_a[] = {
1628 	f2xm1, fyl2x, fptan, fpatan,
1629 	fxtract, fprem1, (FUNC_ST0) fdecstp, (FUNC_ST0) fincstp
1630 };
1631 
1632 void FPU_triga(void)
1633 {
1634 	(trig_table_a[FPU_rm]) (&st(0), FPU_gettag0());
1635 }
1636 
1637 static FUNC_ST0 const trig_table_b[] = {
1638 	fprem, fyl2xp1, fsqrt_, fsincos, frndint_, fscale, (FUNC_ST0) fsin, fcos
1639 };
1640 
1641 void FPU_trigb(void)
1642 {
1643 	(trig_table_b[FPU_rm]) (&st(0), FPU_gettag0());
1644 }
1645