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