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