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