xref: /linux/arch/x86/math-emu/fpu_trig.c (revision 90d32e92011eaae8e70a9169b4e7acf4ca8f9d3a)
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(FPU_REG *st0_ptr, u_char st0_tag)
437 {
438 	clear_C1();
439 	top--;
440 }
441 
442 static void fincstp(FPU_REG *st0_ptr, u_char st0_tag)
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 f_sin(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 void fsin(FPU_REG *st0_ptr, u_char tag)
612 {
613 	f_sin(st0_ptr, tag);
614 }
615 
616 static int f_cos(FPU_REG *st0_ptr, u_char tag)
617 {
618 	u_char st0_sign;
619 
620 	st0_sign = getsign(st0_ptr);
621 
622 	if (tag == TAG_Valid) {
623 		int q;
624 
625 		if (exponent(st0_ptr) > -40) {
626 			if ((exponent(st0_ptr) < 0)
627 			    || ((exponent(st0_ptr) == 0)
628 				&& (significand(st0_ptr) <=
629 				    0xc90fdaa22168c234LL))) {
630 				poly_cos(st0_ptr);
631 
632 				/* We do not really know if up or down */
633 				set_precision_flag_down();
634 
635 				return 0;
636 			} else if ((q = trig_arg(st0_ptr, FCOS)) != -1) {
637 				poly_sine(st0_ptr);
638 
639 				if ((q + 1) & 2)
640 					changesign(st0_ptr);
641 
642 				/* We do not really know if up or down */
643 				set_precision_flag_down();
644 
645 				return 0;
646 			} else {
647 				/* Operand is out of range */
648 				return 1;
649 			}
650 		} else {
651 		      denormal_arg:
652 
653 			setcc(0);
654 			FPU_copy_to_reg0(&CONST_1, TAG_Valid);
655 #ifdef PECULIAR_486
656 			set_precision_flag_down();	/* 80486 appears to do this. */
657 #else
658 			set_precision_flag_up();	/* Must be up. */
659 #endif /* PECULIAR_486 */
660 			return 0;
661 		}
662 	} else if (tag == TAG_Zero) {
663 		FPU_copy_to_reg0(&CONST_1, TAG_Valid);
664 		setcc(0);
665 		return 0;
666 	}
667 
668 	if (tag == TAG_Special)
669 		tag = FPU_Special(st0_ptr);
670 
671 	if (tag == TW_Denormal) {
672 		if (denormal_operand() < 0)
673 			return 1;
674 
675 		goto denormal_arg;
676 	} else if (tag == TW_Infinity) {
677 		/* The 80486 treats infinity as an invalid operand */
678 		arith_invalid(0);
679 		return 1;
680 	} else {
681 		single_arg_error(st0_ptr, tag);	/* requires st0_ptr == &st(0) */
682 		return 1;
683 	}
684 }
685 
686 static void fcos(FPU_REG *st0_ptr, u_char st0_tag)
687 {
688 	f_cos(st0_ptr, st0_tag);
689 }
690 
691 static void fsincos(FPU_REG *st0_ptr, u_char st0_tag)
692 {
693 	FPU_REG *st_new_ptr;
694 	FPU_REG arg;
695 	u_char tag;
696 
697 	/* Stack underflow has higher priority */
698 	if (st0_tag == TAG_Empty) {
699 		FPU_stack_underflow();	/* Puts a QNaN in st(0) */
700 		if (control_word & CW_Invalid) {
701 			st_new_ptr = &st(-1);
702 			push();
703 			FPU_stack_underflow();	/* Puts a QNaN in the new st(0) */
704 		}
705 		return;
706 	}
707 
708 	if (STACK_OVERFLOW) {
709 		FPU_stack_overflow();
710 		return;
711 	}
712 
713 	if (st0_tag == TAG_Special)
714 		tag = FPU_Special(st0_ptr);
715 	else
716 		tag = st0_tag;
717 
718 	if (tag == TW_NaN) {
719 		single_arg_2_error(st0_ptr, TW_NaN);
720 		return;
721 	} else if (tag == TW_Infinity) {
722 		/* The 80486 treats infinity as an invalid operand */
723 		if (arith_invalid(0) >= 0) {
724 			/* Masked response */
725 			push();
726 			arith_invalid(0);
727 		}
728 		return;
729 	}
730 
731 	reg_copy(st0_ptr, &arg);
732 	if (!f_sin(st0_ptr, st0_tag)) {
733 		push();
734 		FPU_copy_to_reg0(&arg, st0_tag);
735 		f_cos(&st(0), st0_tag);
736 	} else {
737 		/* An error, so restore st(0) */
738 		FPU_copy_to_reg0(&arg, st0_tag);
739 	}
740 }
741 
742 /*---------------------------------------------------------------------------*/
743 /* The following all require two arguments: st(0) and st(1) */
744 
745 /* A lean, mean kernel for the fprem instructions. This relies upon
746    the division and rounding to an integer in do_fprem giving an
747    exact result. Because of this, rem_kernel() needs to deal only with
748    the least significant 64 bits, the more significant bits of the
749    result must be zero.
750  */
751 static void rem_kernel(unsigned long long st0, unsigned long long *y,
752 		       unsigned long long st1, unsigned long long q, int n)
753 {
754 	int dummy;
755 	unsigned long long x;
756 
757 	x = st0 << n;
758 
759 	/* Do the required multiplication and subtraction in the one operation */
760 
761 	/* lsw x -= lsw st1 * lsw q */
762 	asm volatile ("mull %4; subl %%eax,%0; sbbl %%edx,%1":"=m"
763 		      (((unsigned *)&x)[0]), "=m"(((unsigned *)&x)[1]),
764 		      "=a"(dummy)
765 		      :"2"(((unsigned *)&st1)[0]), "m"(((unsigned *)&q)[0])
766 		      :"%dx");
767 	/* msw x -= msw st1 * lsw q */
768 	asm volatile ("mull %3; subl %%eax,%0":"=m" (((unsigned *)&x)[1]),
769 		      "=a"(dummy)
770 		      :"1"(((unsigned *)&st1)[1]), "m"(((unsigned *)&q)[0])
771 		      :"%dx");
772 	/* msw x -= lsw st1 * msw q */
773 	asm volatile ("mull %3; subl %%eax,%0":"=m" (((unsigned *)&x)[1]),
774 		      "=a"(dummy)
775 		      :"1"(((unsigned *)&st1)[0]), "m"(((unsigned *)&q)[1])
776 		      :"%dx");
777 
778 	*y = x;
779 }
780 
781 /* Remainder of st(0) / st(1) */
782 /* This routine produces exact results, i.e. there is never any
783    rounding or truncation, etc of the result. */
784 static void do_fprem(FPU_REG *st0_ptr, u_char st0_tag, int round)
785 {
786 	FPU_REG *st1_ptr = &st(1);
787 	u_char st1_tag = FPU_gettagi(1);
788 
789 	if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
790 		FPU_REG tmp, st0, st1;
791 		u_char st0_sign, st1_sign;
792 		u_char tmptag;
793 		int tag;
794 		int old_cw;
795 		int expdif;
796 		long long q;
797 		unsigned short saved_status;
798 		int cc;
799 
800 	      fprem_valid:
801 		/* Convert registers for internal use. */
802 		st0_sign = FPU_to_exp16(st0_ptr, &st0);
803 		st1_sign = FPU_to_exp16(st1_ptr, &st1);
804 		expdif = exponent16(&st0) - exponent16(&st1);
805 
806 		old_cw = control_word;
807 		cc = 0;
808 
809 		/* We want the status following the denorm tests, but don't want
810 		   the status changed by the arithmetic operations. */
811 		saved_status = partial_status;
812 		control_word &= ~CW_RC;
813 		control_word |= RC_CHOP;
814 
815 		if (expdif < 64) {
816 			/* This should be the most common case */
817 
818 			if (expdif > -2) {
819 				u_char sign = st0_sign ^ st1_sign;
820 				tag = FPU_u_div(&st0, &st1, &tmp,
821 						PR_64_BITS | RC_CHOP | 0x3f,
822 						sign);
823 				setsign(&tmp, sign);
824 
825 				if (exponent(&tmp) >= 0) {
826 					FPU_round_to_int(&tmp, tag);	/* Fortunately, this can't
827 									   overflow to 2^64 */
828 					q = significand(&tmp);
829 
830 					rem_kernel(significand(&st0),
831 						   &significand(&tmp),
832 						   significand(&st1),
833 						   q, expdif);
834 
835 					setexponent16(&tmp, exponent16(&st1));
836 				} else {
837 					reg_copy(&st0, &tmp);
838 					q = 0;
839 				}
840 
841 				if ((round == RC_RND)
842 				    && (tmp.sigh & 0xc0000000)) {
843 					/* We may need to subtract st(1) once more,
844 					   to get a result <= 1/2 of st(1). */
845 					unsigned long long x;
846 					expdif =
847 					    exponent16(&st1) - exponent16(&tmp);
848 					if (expdif <= 1) {
849 						if (expdif == 0)
850 							x = significand(&st1) -
851 							    significand(&tmp);
852 						else	/* expdif is 1 */
853 							x = (significand(&st1)
854 							     << 1) -
855 							    significand(&tmp);
856 						if ((x < significand(&tmp)) ||
857 						    /* or equi-distant (from 0 & st(1)) and q is odd */
858 						    ((x == significand(&tmp))
859 						     && (q & 1))) {
860 							st0_sign = !st0_sign;
861 							significand(&tmp) = x;
862 							q++;
863 						}
864 					}
865 				}
866 
867 				if (q & 4)
868 					cc |= SW_C0;
869 				if (q & 2)
870 					cc |= SW_C3;
871 				if (q & 1)
872 					cc |= SW_C1;
873 			} else {
874 				control_word = old_cw;
875 				setcc(0);
876 				return;
877 			}
878 		} else {
879 			/* There is a large exponent difference ( >= 64 ) */
880 			/* To make much sense, the code in this section should
881 			   be done at high precision. */
882 			int exp_1, N;
883 			u_char sign;
884 
885 			/* prevent overflow here */
886 			/* N is 'a number between 32 and 63' (p26-113) */
887 			reg_copy(&st0, &tmp);
888 			tmptag = st0_tag;
889 			N = (expdif & 0x0000001f) + 32;	/* This choice gives results
890 							   identical to an AMD 486 */
891 			setexponent16(&tmp, N);
892 			exp_1 = exponent16(&st1);
893 			setexponent16(&st1, 0);
894 			expdif -= N;
895 
896 			sign = getsign(&tmp) ^ st1_sign;
897 			tag =
898 			    FPU_u_div(&tmp, &st1, &tmp,
899 				      PR_64_BITS | RC_CHOP | 0x3f, sign);
900 			setsign(&tmp, sign);
901 
902 			FPU_round_to_int(&tmp, tag);	/* Fortunately, this can't
903 							   overflow to 2^64 */
904 
905 			rem_kernel(significand(&st0),
906 				   &significand(&tmp),
907 				   significand(&st1),
908 				   significand(&tmp), exponent(&tmp)
909 			    );
910 			setexponent16(&tmp, exp_1 + expdif);
911 
912 			/* It is possible for the operation to be complete here.
913 			   What does the IEEE standard say? The Intel 80486 manual
914 			   implies that the operation will never be completed at this
915 			   point, and the behaviour of a real 80486 confirms this.
916 			 */
917 			if (!(tmp.sigh | tmp.sigl)) {
918 				/* The result is zero */
919 				control_word = old_cw;
920 				partial_status = saved_status;
921 				FPU_copy_to_reg0(&CONST_Z, TAG_Zero);
922 				setsign(&st0, st0_sign);
923 #ifdef PECULIAR_486
924 				setcc(SW_C2);
925 #else
926 				setcc(0);
927 #endif /* PECULIAR_486 */
928 				return;
929 			}
930 			cc = SW_C2;
931 		}
932 
933 		control_word = old_cw;
934 		partial_status = saved_status;
935 		tag = FPU_normalize_nuo(&tmp);
936 		reg_copy(&tmp, st0_ptr);
937 
938 		/* The only condition to be looked for is underflow,
939 		   and it can occur here only if underflow is unmasked. */
940 		if ((exponent16(&tmp) <= EXP_UNDER) && (tag != TAG_Zero)
941 		    && !(control_word & CW_Underflow)) {
942 			setcc(cc);
943 			tag = arith_underflow(st0_ptr);
944 			setsign(st0_ptr, st0_sign);
945 			FPU_settag0(tag);
946 			return;
947 		} else if ((exponent16(&tmp) > EXP_UNDER) || (tag == TAG_Zero)) {
948 			stdexp(st0_ptr);
949 			setsign(st0_ptr, st0_sign);
950 		} else {
951 			tag =
952 			    FPU_round(st0_ptr, 0, 0, FULL_PRECISION, st0_sign);
953 		}
954 		FPU_settag0(tag);
955 		setcc(cc);
956 
957 		return;
958 	}
959 
960 	if (st0_tag == TAG_Special)
961 		st0_tag = FPU_Special(st0_ptr);
962 	if (st1_tag == TAG_Special)
963 		st1_tag = FPU_Special(st1_ptr);
964 
965 	if (((st0_tag == TAG_Valid) && (st1_tag == TW_Denormal))
966 	    || ((st0_tag == TW_Denormal) && (st1_tag == TAG_Valid))
967 	    || ((st0_tag == TW_Denormal) && (st1_tag == TW_Denormal))) {
968 		if (denormal_operand() < 0)
969 			return;
970 		goto fprem_valid;
971 	} else if ((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty)) {
972 		FPU_stack_underflow();
973 		return;
974 	} else if (st0_tag == TAG_Zero) {
975 		if (st1_tag == TAG_Valid) {
976 			setcc(0);
977 			return;
978 		} else if (st1_tag == TW_Denormal) {
979 			if (denormal_operand() < 0)
980 				return;
981 			setcc(0);
982 			return;
983 		} else if (st1_tag == TAG_Zero) {
984 			arith_invalid(0);
985 			return;
986 		} /* fprem(?,0) always invalid */
987 		else if (st1_tag == TW_Infinity) {
988 			setcc(0);
989 			return;
990 		}
991 	} else if ((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal)) {
992 		if (st1_tag == TAG_Zero) {
993 			arith_invalid(0);	/* fprem(Valid,Zero) is invalid */
994 			return;
995 		} else if (st1_tag != TW_NaN) {
996 			if (((st0_tag == TW_Denormal)
997 			     || (st1_tag == TW_Denormal))
998 			    && (denormal_operand() < 0))
999 				return;
1000 
1001 			if (st1_tag == TW_Infinity) {
1002 				/* fprem(Valid,Infinity) is o.k. */
1003 				setcc(0);
1004 				return;
1005 			}
1006 		}
1007 	} else if (st0_tag == TW_Infinity) {
1008 		if (st1_tag != TW_NaN) {
1009 			arith_invalid(0);	/* fprem(Infinity,?) is invalid */
1010 			return;
1011 		}
1012 	}
1013 
1014 	/* One of the registers must contain a NaN if we got here. */
1015 
1016 #ifdef PARANOID
1017 	if ((st0_tag != TW_NaN) && (st1_tag != TW_NaN))
1018 		EXCEPTION(EX_INTERNAL | 0x118);
1019 #endif /* PARANOID */
1020 
1021 	real_2op_NaN(st1_ptr, st1_tag, 0, st1_ptr);
1022 
1023 }
1024 
1025 /* ST(1) <- ST(1) * log ST;  pop ST */
1026 static void fyl2x(FPU_REG *st0_ptr, u_char st0_tag)
1027 {
1028 	FPU_REG *st1_ptr = &st(1), exponent;
1029 	u_char st1_tag = FPU_gettagi(1);
1030 	u_char sign;
1031 	int e, tag;
1032 
1033 	clear_C1();
1034 
1035 	if ((st0_tag == TAG_Valid) && (st1_tag == TAG_Valid)) {
1036 	      both_valid:
1037 		/* Both regs are Valid or Denormal */
1038 		if (signpositive(st0_ptr)) {
1039 			if (st0_tag == TW_Denormal)
1040 				FPU_to_exp16(st0_ptr, st0_ptr);
1041 			else
1042 				/* Convert st(0) for internal use. */
1043 				setexponent16(st0_ptr, exponent(st0_ptr));
1044 
1045 			if ((st0_ptr->sigh == 0x80000000)
1046 			    && (st0_ptr->sigl == 0)) {
1047 				/* Special case. The result can be precise. */
1048 				u_char esign;
1049 				e = exponent16(st0_ptr);
1050 				if (e >= 0) {
1051 					exponent.sigh = e;
1052 					esign = SIGN_POS;
1053 				} else {
1054 					exponent.sigh = -e;
1055 					esign = SIGN_NEG;
1056 				}
1057 				exponent.sigl = 0;
1058 				setexponent16(&exponent, 31);
1059 				tag = FPU_normalize_nuo(&exponent);
1060 				stdexp(&exponent);
1061 				setsign(&exponent, esign);
1062 				tag =
1063 				    FPU_mul(&exponent, tag, 1, FULL_PRECISION);
1064 				if (tag >= 0)
1065 					FPU_settagi(1, tag);
1066 			} else {
1067 				/* The usual case */
1068 				sign = getsign(st1_ptr);
1069 				if (st1_tag == TW_Denormal)
1070 					FPU_to_exp16(st1_ptr, st1_ptr);
1071 				else
1072 					/* Convert st(1) for internal use. */
1073 					setexponent16(st1_ptr,
1074 						      exponent(st1_ptr));
1075 				poly_l2(st0_ptr, st1_ptr, sign);
1076 			}
1077 		} else {
1078 			/* negative */
1079 			if (arith_invalid(1) < 0)
1080 				return;
1081 		}
1082 
1083 		FPU_pop();
1084 
1085 		return;
1086 	}
1087 
1088 	if (st0_tag == TAG_Special)
1089 		st0_tag = FPU_Special(st0_ptr);
1090 	if (st1_tag == TAG_Special)
1091 		st1_tag = FPU_Special(st1_ptr);
1092 
1093 	if ((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty)) {
1094 		FPU_stack_underflow_pop(1);
1095 		return;
1096 	} else if ((st0_tag <= TW_Denormal) && (st1_tag <= TW_Denormal)) {
1097 		if (st0_tag == TAG_Zero) {
1098 			if (st1_tag == TAG_Zero) {
1099 				/* Both args zero is invalid */
1100 				if (arith_invalid(1) < 0)
1101 					return;
1102 			} else {
1103 				u_char sign;
1104 				sign = getsign(st1_ptr) ^ SIGN_NEG;
1105 				if (FPU_divide_by_zero(1, sign) < 0)
1106 					return;
1107 
1108 				setsign(st1_ptr, sign);
1109 			}
1110 		} else if (st1_tag == TAG_Zero) {
1111 			/* st(1) contains zero, st(0) valid <> 0 */
1112 			/* Zero is the valid answer */
1113 			sign = getsign(st1_ptr);
1114 
1115 			if (signnegative(st0_ptr)) {
1116 				/* log(negative) */
1117 				if (arith_invalid(1) < 0)
1118 					return;
1119 			} else if ((st0_tag == TW_Denormal)
1120 				   && (denormal_operand() < 0))
1121 				return;
1122 			else {
1123 				if (exponent(st0_ptr) < 0)
1124 					sign ^= SIGN_NEG;
1125 
1126 				FPU_copy_to_reg1(&CONST_Z, TAG_Zero);
1127 				setsign(st1_ptr, sign);
1128 			}
1129 		} else {
1130 			/* One or both operands are denormals. */
1131 			if (denormal_operand() < 0)
1132 				return;
1133 			goto both_valid;
1134 		}
1135 	} else if ((st0_tag == TW_NaN) || (st1_tag == TW_NaN)) {
1136 		if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1137 			return;
1138 	}
1139 	/* One or both arg must be an infinity */
1140 	else if (st0_tag == TW_Infinity) {
1141 		if ((signnegative(st0_ptr)) || (st1_tag == TAG_Zero)) {
1142 			/* log(-infinity) or 0*log(infinity) */
1143 			if (arith_invalid(1) < 0)
1144 				return;
1145 		} else {
1146 			u_char sign = getsign(st1_ptr);
1147 
1148 			if ((st1_tag == TW_Denormal)
1149 			    && (denormal_operand() < 0))
1150 				return;
1151 
1152 			FPU_copy_to_reg1(&CONST_INF, TAG_Special);
1153 			setsign(st1_ptr, sign);
1154 		}
1155 	}
1156 	/* st(1) must be infinity here */
1157 	else if (((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal))
1158 		 && (signpositive(st0_ptr))) {
1159 		if (exponent(st0_ptr) >= 0) {
1160 			if ((exponent(st0_ptr) == 0) &&
1161 			    (st0_ptr->sigh == 0x80000000) &&
1162 			    (st0_ptr->sigl == 0)) {
1163 				/* st(0) holds 1.0 */
1164 				/* infinity*log(1) */
1165 				if (arith_invalid(1) < 0)
1166 					return;
1167 			}
1168 			/* else st(0) is positive and > 1.0 */
1169 		} else {
1170 			/* st(0) is positive and < 1.0 */
1171 
1172 			if ((st0_tag == TW_Denormal)
1173 			    && (denormal_operand() < 0))
1174 				return;
1175 
1176 			changesign(st1_ptr);
1177 		}
1178 	} else {
1179 		/* st(0) must be zero or negative */
1180 		if (st0_tag == TAG_Zero) {
1181 			/* This should be invalid, but a real 80486 is happy with it. */
1182 
1183 #ifndef PECULIAR_486
1184 			sign = getsign(st1_ptr);
1185 			if (FPU_divide_by_zero(1, sign) < 0)
1186 				return;
1187 #endif /* PECULIAR_486 */
1188 
1189 			changesign(st1_ptr);
1190 		} else if (arith_invalid(1) < 0)	/* log(negative) */
1191 			return;
1192 	}
1193 
1194 	FPU_pop();
1195 }
1196 
1197 static void fpatan(FPU_REG *st0_ptr, u_char st0_tag)
1198 {
1199 	FPU_REG *st1_ptr = &st(1);
1200 	u_char st1_tag = FPU_gettagi(1);
1201 	int tag;
1202 
1203 	clear_C1();
1204 	if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
1205 	      valid_atan:
1206 
1207 		poly_atan(st0_ptr, st0_tag, st1_ptr, st1_tag);
1208 
1209 		FPU_pop();
1210 
1211 		return;
1212 	}
1213 
1214 	if (st0_tag == TAG_Special)
1215 		st0_tag = FPU_Special(st0_ptr);
1216 	if (st1_tag == TAG_Special)
1217 		st1_tag = FPU_Special(st1_ptr);
1218 
1219 	if (((st0_tag == TAG_Valid) && (st1_tag == TW_Denormal))
1220 	    || ((st0_tag == TW_Denormal) && (st1_tag == TAG_Valid))
1221 	    || ((st0_tag == TW_Denormal) && (st1_tag == TW_Denormal))) {
1222 		if (denormal_operand() < 0)
1223 			return;
1224 
1225 		goto valid_atan;
1226 	} else if ((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty)) {
1227 		FPU_stack_underflow_pop(1);
1228 		return;
1229 	} else if ((st0_tag == TW_NaN) || (st1_tag == TW_NaN)) {
1230 		if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) >= 0)
1231 			FPU_pop();
1232 		return;
1233 	} else if ((st0_tag == TW_Infinity) || (st1_tag == TW_Infinity)) {
1234 		u_char sign = getsign(st1_ptr);
1235 		if (st0_tag == TW_Infinity) {
1236 			if (st1_tag == TW_Infinity) {
1237 				if (signpositive(st0_ptr)) {
1238 					FPU_copy_to_reg1(&CONST_PI4, TAG_Valid);
1239 				} else {
1240 					setpositive(st1_ptr);
1241 					tag =
1242 					    FPU_u_add(&CONST_PI4, &CONST_PI2,
1243 						      st1_ptr, FULL_PRECISION,
1244 						      SIGN_POS,
1245 						      exponent(&CONST_PI4),
1246 						      exponent(&CONST_PI2));
1247 					if (tag >= 0)
1248 						FPU_settagi(1, tag);
1249 				}
1250 			} else {
1251 				if ((st1_tag == TW_Denormal)
1252 				    && (denormal_operand() < 0))
1253 					return;
1254 
1255 				if (signpositive(st0_ptr)) {
1256 					FPU_copy_to_reg1(&CONST_Z, TAG_Zero);
1257 					setsign(st1_ptr, sign);	/* An 80486 preserves the sign */
1258 					FPU_pop();
1259 					return;
1260 				} else {
1261 					FPU_copy_to_reg1(&CONST_PI, TAG_Valid);
1262 				}
1263 			}
1264 		} else {
1265 			/* st(1) is infinity, st(0) not infinity */
1266 			if ((st0_tag == TW_Denormal)
1267 			    && (denormal_operand() < 0))
1268 				return;
1269 
1270 			FPU_copy_to_reg1(&CONST_PI2, TAG_Valid);
1271 		}
1272 		setsign(st1_ptr, sign);
1273 	} else if (st1_tag == TAG_Zero) {
1274 		/* st(0) must be valid or zero */
1275 		u_char sign = getsign(st1_ptr);
1276 
1277 		if ((st0_tag == TW_Denormal) && (denormal_operand() < 0))
1278 			return;
1279 
1280 		if (signpositive(st0_ptr)) {
1281 			/* An 80486 preserves the sign */
1282 			FPU_pop();
1283 			return;
1284 		}
1285 
1286 		FPU_copy_to_reg1(&CONST_PI, TAG_Valid);
1287 		setsign(st1_ptr, sign);
1288 	} else if (st0_tag == TAG_Zero) {
1289 		/* st(1) must be TAG_Valid here */
1290 		u_char sign = getsign(st1_ptr);
1291 
1292 		if ((st1_tag == TW_Denormal) && (denormal_operand() < 0))
1293 			return;
1294 
1295 		FPU_copy_to_reg1(&CONST_PI2, TAG_Valid);
1296 		setsign(st1_ptr, sign);
1297 	}
1298 #ifdef PARANOID
1299 	else
1300 		EXCEPTION(EX_INTERNAL | 0x125);
1301 #endif /* PARANOID */
1302 
1303 	FPU_pop();
1304 	set_precision_flag_up();	/* We do not really know if up or down */
1305 }
1306 
1307 static void fprem(FPU_REG *st0_ptr, u_char st0_tag)
1308 {
1309 	do_fprem(st0_ptr, st0_tag, RC_CHOP);
1310 }
1311 
1312 static void fprem1(FPU_REG *st0_ptr, u_char st0_tag)
1313 {
1314 	do_fprem(st0_ptr, st0_tag, RC_RND);
1315 }
1316 
1317 static void fyl2xp1(FPU_REG *st0_ptr, u_char st0_tag)
1318 {
1319 	u_char sign, sign1;
1320 	FPU_REG *st1_ptr = &st(1), a, b;
1321 	u_char st1_tag = FPU_gettagi(1);
1322 
1323 	clear_C1();
1324 	if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
1325 	      valid_yl2xp1:
1326 
1327 		sign = getsign(st0_ptr);
1328 		sign1 = getsign(st1_ptr);
1329 
1330 		FPU_to_exp16(st0_ptr, &a);
1331 		FPU_to_exp16(st1_ptr, &b);
1332 
1333 		if (poly_l2p1(sign, sign1, &a, &b, st1_ptr))
1334 			return;
1335 
1336 		FPU_pop();
1337 		return;
1338 	}
1339 
1340 	if (st0_tag == TAG_Special)
1341 		st0_tag = FPU_Special(st0_ptr);
1342 	if (st1_tag == TAG_Special)
1343 		st1_tag = FPU_Special(st1_ptr);
1344 
1345 	if (((st0_tag == TAG_Valid) && (st1_tag == TW_Denormal))
1346 	    || ((st0_tag == TW_Denormal) && (st1_tag == TAG_Valid))
1347 	    || ((st0_tag == TW_Denormal) && (st1_tag == TW_Denormal))) {
1348 		if (denormal_operand() < 0)
1349 			return;
1350 
1351 		goto valid_yl2xp1;
1352 	} else if ((st0_tag == TAG_Empty) | (st1_tag == TAG_Empty)) {
1353 		FPU_stack_underflow_pop(1);
1354 		return;
1355 	} else if (st0_tag == TAG_Zero) {
1356 		switch (st1_tag) {
1357 		case TW_Denormal:
1358 			if (denormal_operand() < 0)
1359 				return;
1360 			fallthrough;
1361 		case TAG_Zero:
1362 		case TAG_Valid:
1363 			setsign(st0_ptr, getsign(st0_ptr) ^ getsign(st1_ptr));
1364 			FPU_copy_to_reg1(st0_ptr, st0_tag);
1365 			break;
1366 
1367 		case TW_Infinity:
1368 			/* Infinity*log(1) */
1369 			if (arith_invalid(1) < 0)
1370 				return;
1371 			break;
1372 
1373 		case TW_NaN:
1374 			if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1375 				return;
1376 			break;
1377 
1378 		default:
1379 #ifdef PARANOID
1380 			EXCEPTION(EX_INTERNAL | 0x116);
1381 			return;
1382 #endif /* PARANOID */
1383 			break;
1384 		}
1385 	} else if ((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal)) {
1386 		switch (st1_tag) {
1387 		case TAG_Zero:
1388 			if (signnegative(st0_ptr)) {
1389 				if (exponent(st0_ptr) >= 0) {
1390 					/* st(0) holds <= -1.0 */
1391 #ifdef PECULIAR_486		/* Stupid 80486 doesn't worry about log(negative). */
1392 					changesign(st1_ptr);
1393 #else
1394 					if (arith_invalid(1) < 0)
1395 						return;
1396 #endif /* PECULIAR_486 */
1397 				} else if ((st0_tag == TW_Denormal)
1398 					   && (denormal_operand() < 0))
1399 					return;
1400 				else
1401 					changesign(st1_ptr);
1402 			} else if ((st0_tag == TW_Denormal)
1403 				   && (denormal_operand() < 0))
1404 				return;
1405 			break;
1406 
1407 		case TW_Infinity:
1408 			if (signnegative(st0_ptr)) {
1409 				if ((exponent(st0_ptr) >= 0) &&
1410 				    !((st0_ptr->sigh == 0x80000000) &&
1411 				      (st0_ptr->sigl == 0))) {
1412 					/* st(0) holds < -1.0 */
1413 #ifdef PECULIAR_486		/* Stupid 80486 doesn't worry about log(negative). */
1414 					changesign(st1_ptr);
1415 #else
1416 					if (arith_invalid(1) < 0)
1417 						return;
1418 #endif /* PECULIAR_486 */
1419 				} else if ((st0_tag == TW_Denormal)
1420 					   && (denormal_operand() < 0))
1421 					return;
1422 				else
1423 					changesign(st1_ptr);
1424 			} else if ((st0_tag == TW_Denormal)
1425 				   && (denormal_operand() < 0))
1426 				return;
1427 			break;
1428 
1429 		case TW_NaN:
1430 			if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1431 				return;
1432 		}
1433 
1434 	} else if (st0_tag == TW_NaN) {
1435 		if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1436 			return;
1437 	} else if (st0_tag == TW_Infinity) {
1438 		if (st1_tag == TW_NaN) {
1439 			if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1440 				return;
1441 		} else if (signnegative(st0_ptr)) {
1442 #ifndef PECULIAR_486
1443 			/* This should have higher priority than denormals, but... */
1444 			if (arith_invalid(1) < 0)	/* log(-infinity) */
1445 				return;
1446 #endif /* PECULIAR_486 */
1447 			if ((st1_tag == TW_Denormal)
1448 			    && (denormal_operand() < 0))
1449 				return;
1450 #ifdef PECULIAR_486
1451 			/* Denormal operands actually get higher priority */
1452 			if (arith_invalid(1) < 0)	/* log(-infinity) */
1453 				return;
1454 #endif /* PECULIAR_486 */
1455 		} else if (st1_tag == TAG_Zero) {
1456 			/* log(infinity) */
1457 			if (arith_invalid(1) < 0)
1458 				return;
1459 		}
1460 
1461 		/* st(1) must be valid here. */
1462 
1463 		else if ((st1_tag == TW_Denormal) && (denormal_operand() < 0))
1464 			return;
1465 
1466 		/* The Manual says that log(Infinity) is invalid, but a real
1467 		   80486 sensibly says that it is o.k. */
1468 		else {
1469 			u_char sign = getsign(st1_ptr);
1470 			FPU_copy_to_reg1(&CONST_INF, TAG_Special);
1471 			setsign(st1_ptr, sign);
1472 		}
1473 	}
1474 #ifdef PARANOID
1475 	else {
1476 		EXCEPTION(EX_INTERNAL | 0x117);
1477 		return;
1478 	}
1479 #endif /* PARANOID */
1480 
1481 	FPU_pop();
1482 	return;
1483 
1484 }
1485 
1486 static void fscale(FPU_REG *st0_ptr, u_char st0_tag)
1487 {
1488 	FPU_REG *st1_ptr = &st(1);
1489 	u_char st1_tag = FPU_gettagi(1);
1490 	int old_cw = control_word;
1491 	u_char sign = getsign(st0_ptr);
1492 
1493 	clear_C1();
1494 	if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
1495 		long scale;
1496 		FPU_REG tmp;
1497 
1498 		/* Convert register for internal use. */
1499 		setexponent16(st0_ptr, exponent(st0_ptr));
1500 
1501 	      valid_scale:
1502 
1503 		if (exponent(st1_ptr) > 30) {
1504 			/* 2^31 is far too large, would require 2^(2^30) or 2^(-2^30) */
1505 
1506 			if (signpositive(st1_ptr)) {
1507 				EXCEPTION(EX_Overflow);
1508 				FPU_copy_to_reg0(&CONST_INF, TAG_Special);
1509 			} else {
1510 				EXCEPTION(EX_Underflow);
1511 				FPU_copy_to_reg0(&CONST_Z, TAG_Zero);
1512 			}
1513 			setsign(st0_ptr, sign);
1514 			return;
1515 		}
1516 
1517 		control_word &= ~CW_RC;
1518 		control_word |= RC_CHOP;
1519 		reg_copy(st1_ptr, &tmp);
1520 		FPU_round_to_int(&tmp, st1_tag);	/* This can never overflow here */
1521 		control_word = old_cw;
1522 		scale = signnegative(st1_ptr) ? -tmp.sigl : tmp.sigl;
1523 		scale += exponent16(st0_ptr);
1524 
1525 		setexponent16(st0_ptr, scale);
1526 
1527 		/* Use FPU_round() to properly detect under/overflow etc */
1528 		FPU_round(st0_ptr, 0, 0, control_word, sign);
1529 
1530 		return;
1531 	}
1532 
1533 	if (st0_tag == TAG_Special)
1534 		st0_tag = FPU_Special(st0_ptr);
1535 	if (st1_tag == TAG_Special)
1536 		st1_tag = FPU_Special(st1_ptr);
1537 
1538 	if ((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal)) {
1539 		switch (st1_tag) {
1540 		case TAG_Valid:
1541 			/* st(0) must be a denormal */
1542 			if ((st0_tag == TW_Denormal)
1543 			    && (denormal_operand() < 0))
1544 				return;
1545 
1546 			FPU_to_exp16(st0_ptr, st0_ptr);	/* Will not be left on stack */
1547 			goto valid_scale;
1548 
1549 		case TAG_Zero:
1550 			if (st0_tag == TW_Denormal)
1551 				denormal_operand();
1552 			return;
1553 
1554 		case TW_Denormal:
1555 			denormal_operand();
1556 			return;
1557 
1558 		case TW_Infinity:
1559 			if ((st0_tag == TW_Denormal)
1560 			    && (denormal_operand() < 0))
1561 				return;
1562 
1563 			if (signpositive(st1_ptr))
1564 				FPU_copy_to_reg0(&CONST_INF, TAG_Special);
1565 			else
1566 				FPU_copy_to_reg0(&CONST_Z, TAG_Zero);
1567 			setsign(st0_ptr, sign);
1568 			return;
1569 
1570 		case TW_NaN:
1571 			real_2op_NaN(st1_ptr, st1_tag, 0, st0_ptr);
1572 			return;
1573 		}
1574 	} else if (st0_tag == TAG_Zero) {
1575 		switch (st1_tag) {
1576 		case TAG_Valid:
1577 		case TAG_Zero:
1578 			return;
1579 
1580 		case TW_Denormal:
1581 			denormal_operand();
1582 			return;
1583 
1584 		case TW_Infinity:
1585 			if (signpositive(st1_ptr))
1586 				arith_invalid(0);	/* Zero scaled by +Infinity */
1587 			return;
1588 
1589 		case TW_NaN:
1590 			real_2op_NaN(st1_ptr, st1_tag, 0, st0_ptr);
1591 			return;
1592 		}
1593 	} else if (st0_tag == TW_Infinity) {
1594 		switch (st1_tag) {
1595 		case TAG_Valid:
1596 		case TAG_Zero:
1597 			return;
1598 
1599 		case TW_Denormal:
1600 			denormal_operand();
1601 			return;
1602 
1603 		case TW_Infinity:
1604 			if (signnegative(st1_ptr))
1605 				arith_invalid(0);	/* Infinity scaled by -Infinity */
1606 			return;
1607 
1608 		case TW_NaN:
1609 			real_2op_NaN(st1_ptr, st1_tag, 0, st0_ptr);
1610 			return;
1611 		}
1612 	} else if (st0_tag == TW_NaN) {
1613 		if (st1_tag != TAG_Empty) {
1614 			real_2op_NaN(st1_ptr, st1_tag, 0, st0_ptr);
1615 			return;
1616 		}
1617 	}
1618 #ifdef PARANOID
1619 	if (!((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty))) {
1620 		EXCEPTION(EX_INTERNAL | 0x115);
1621 		return;
1622 	}
1623 #endif
1624 
1625 	/* At least one of st(0), st(1) must be empty */
1626 	FPU_stack_underflow();
1627 
1628 }
1629 
1630 /*---------------------------------------------------------------------------*/
1631 
1632 static FUNC_ST0 const trig_table_a[] = {
1633 	f2xm1, fyl2x, fptan, fpatan,
1634 	fxtract, fprem1, fdecstp, fincstp,
1635 };
1636 
1637 void FPU_triga(void)
1638 {
1639 	(trig_table_a[FPU_rm]) (&st(0), FPU_gettag0());
1640 }
1641 
1642 static FUNC_ST0 const trig_table_b[] = {
1643 	fprem, fyl2xp1, fsqrt_, fsincos, frndint_, fscale, fsin, fcos
1644 };
1645 
1646 void FPU_trigb(void)
1647 {
1648 	(trig_table_b[FPU_rm]) (&st(0), FPU_gettag0());
1649 }
1650