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 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
611static void fsin(FPU_REG *st0_ptr, u_char tag)
612{
613	f_sin(st0_ptr, tag);
614}
615
616static 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
686static void fcos(FPU_REG *st0_ptr, u_char st0_tag)
687{
688	f_cos(st0_ptr, st0_tag);
689}
690
691static 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 */
751static 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. */
784static 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 */
1026static 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
1197static 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
1307static void fprem(FPU_REG *st0_ptr, u_char st0_tag)
1308{
1309	do_fprem(st0_ptr, st0_tag, RC_CHOP);
1310}
1311
1312static void fprem1(FPU_REG *st0_ptr, u_char st0_tag)
1313{
1314	do_fprem(st0_ptr, st0_tag, RC_RND);
1315}
1316
1317static 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
1486static 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
1632static FUNC_ST0 const trig_table_a[] = {
1633	f2xm1, fyl2x, fptan, fpatan,
1634	fxtract, fprem1, (FUNC_ST0) fdecstp, (FUNC_ST0) fincstp
1635};
1636
1637void FPU_triga(void)
1638{
1639	(trig_table_a[FPU_rm]) (&st(0), FPU_gettag0());
1640}
1641
1642static FUNC_ST0 const trig_table_b[] = {
1643	fprem, fyl2xp1, fsqrt_, fsincos, frndint_, fscale, fsin, fcos
1644};
1645
1646void FPU_trigb(void)
1647{
1648	(trig_table_b[FPU_rm]) (&st(0), FPU_gettag0());
1649}
1650