GCC Wikia
Advertisement

このページを編集する際は,編集に関する方針に従ってください.[]

概要[]

引数[]

実装[]

555 /* Divide doubleword integer LNUM, HNUM by doubleword integer LDEN, HDEN
556    for a quotient (stored in *LQUO, *HQUO) and remainder (in *LREM, *HREM).
557    CODE is a tree code for a kind of division, one of
558    TRUNC_DIV_EXPR, FLOOR_DIV_EXPR, CEIL_DIV_EXPR, ROUND_DIV_EXPR
559    or EXACT_DIV_EXPR
560    It controls how the quotient is rounded to an integer.
561    Return nonzero if the operation overflows.
562    UNS nonzero says do unsigned division.  */
563 
564 int
565 div_and_round_double (enum tree_code code, int uns,
566                       unsigned HOST_WIDE_INT lnum_orig, /* num == numerator == dividend */
567                       HOST_WIDE_INT hnum_orig,
568                       unsigned HOST_WIDE_INT lden_orig, /* den == denominator == divisor */
569                       HOST_WIDE_INT hden_orig,
570                       unsigned HOST_WIDE_INT *lquo,
571                       HOST_WIDE_INT *hquo, unsigned HOST_WIDE_INT *lrem,
572                       HOST_WIDE_INT *hrem)
573 {
574   int quo_neg = 0;
575   HOST_WIDE_INT num[4 + 1];     /* extra element for scaling.  */
576   HOST_WIDE_INT den[4], quo[4];
577   int i, j;
578   unsigned HOST_WIDE_INT work;
579   unsigned HOST_WIDE_INT carry = 0;
580   unsigned HOST_WIDE_INT lnum = lnum_orig;
581   HOST_WIDE_INT hnum = hnum_orig;
582   unsigned HOST_WIDE_INT lden = lden_orig;
583   HOST_WIDE_INT hden = hden_orig;
584   int overflow = 0;
585 

  • ゼロ除算

586   if (hden == 0 && lden == 0)
587     overflow = 1, lden = 1;
588 

  • 符号の計算

589   /* Calculate quotient sign and convert operands to unsigned.  */
590   if (!uns)
591     {
592       if (hnum < 0)
593         {
594           quo_neg = ~ quo_neg;
595           /* (minimum integer) / (-1) is the only overflow case.  */
596           if (neg_double (lnum, hnum, &lnum, &hnum)
597               && ((HOST_WIDE_INT) lden & hden) == -1)
598             overflow = 1;
599         }
600       if (hden < 0)
601         {
602           quo_neg = ~ quo_neg;
603           neg_double (lden, hden, &lden, &hden);
604         }
605     }
606 

  • 上部がどちらもゼロ
    • 下部同士の計算で済む

607   if (hnum == 0 && hden == 0)
608     {                           /* single precision */
609       *hquo = *hrem = 0;
610       /* This unsigned division rounds toward zero.  */
611       *lquo = lnum / lden;
612       goto finish_up;
613     }
614 

  • 割られる数の方が小さい

615   if (hnum == 0)
616     {                           /* trivial case: dividend < divisor */
617       /* hden != 0 already checked.  */
618       *hquo = *lquo = 0;
619       *hrem = hnum;
620       *lrem = lnum;
621       goto finish_up;
622     }
623 


  • 初期化

624   memset (quo, 0, sizeof quo);
625 
626   memset (num, 0, sizeof num);  /* to zero 9th element */
627   memset (den, 0, sizeof den);
628 
629   encode (num, lnum, hnum);
630   encode (den, lden, hden);
631 

  • ldenは最大値の半分以下
  • 筆算方式
    • q0がquo[0], n0がnum[0]
      q3 q2 q1 q0
      ___________
lden /n3 n2 n1 n0
      XX
      -----------
      carry

632   /* Special code for when the divisor < BASE.  */
633   if (hden == 0 && lden < (unsigned HOST_WIDE_INT) BASE)
634     {
635       /* hnum != 0 already checked.  */
636       for (i = 4 - 1; i >= 0; i--)
637         {
638           work = num[i] + carry * BASE;
639           quo[i] = work / lden;
640           carry = work % lden;
641         }
642     }


643   else
644     {
645       /* Full double precision division,
646          with thanks to Don Knuth's "Seminumerical Algorithms".  */
647       int num_hi_sig, den_hi_sig;
648       unsigned HOST_WIDE_INT quo_est, scale;
649

  • どこまで割る数か調べる

650       /* Find the highest nonzero divisor digit.  */
651       for (i = 4 - 1;; i--)
652         if (den[i] != 0)
653           {
654             den_hi_sig = i;
655             break;
656           }
657 
658       /* Insure that the first digit of the divisor is at least BASE/2.
659          This is required by the quotient digit estimation algorithm.  */
660 

  • den[den_hi_sig]が負の数でなければ1より大きくなる
  • 割られる数(num), 割る数(den)をscale倍する
    • den[den_hi_sig]の最上位ビットを1にする必要がある?

661       scale = BASE / (den[den_hi_sig] + 1);
662       if (scale > 1)
663         {               /* scale divisor and dividend */
664           carry = 0;
665           for (i = 0; i <= 4 - 1; i++)
666             {
667               work = (num[i] * scale) + carry;
668               num[i] = LOWPART (work);
669               carry = HIGHPART (work);
670             }
671 
672           num[4] = carry;
673           carry = 0;
674           for (i = 0; i <= 4 - 1; i++)
675             {
676               work = (den[i] * scale) + carry;
677               den[i] = LOWPART (work);
678               carry = HIGHPART (work);
679               if (den[i] != 0) den_hi_sig = i;
680             }
681         }
682 


683       num_hi_sig = 4;
684 
685       /* Main loop */
686       for (i = num_hi_sig - den_hi_sig - 1; i >= 0; i--)
687         {
688           /* Guess the next quotient digit, quo_est, by dividing the first
689              two remaining dividend digits by the high order quotient digit.
690              quo_est is never low and is at most 2 high.  */
691           unsigned HOST_WIDE_INT tmp;
692 
693           num_hi_sig = i + den_hi_sig + 1;
694           work = num[num_hi_sig] * BASE + num[num_hi_sig - 1];
695           if (num[num_hi_sig] != den[den_hi_sig])
696             quo_est = work / den[den_hi_sig];
697           else
698             quo_est = BASE - 1;
699 
700           /* Refine quo_est so it's usually correct, and at most one high.  */
701           tmp = work - quo_est * den[den_hi_sig];
702           if (tmp < BASE
703               && (den[den_hi_sig - 1] * quo_est
704                   > (tmp * BASE + num[num_hi_sig - 2])))
705             quo_est--;
706 
707           /* Try QUO_EST as the quotient digit, by multiplying the
708              divisor by QUO_EST and subtracting from the remaining dividend.
709              Keep in mind that QUO_EST is the I - 1st digit.  */
710 
711           carry = 0;
712           for (j = 0; j <= den_hi_sig; j++)
713             {
714               work = quo_est * den[j] + carry;
715               carry = HIGHPART (work);
716               work = num[i + j] - LOWPART (work);
717               num[i + j] = LOWPART (work);
718               carry += HIGHPART (work) != 0;
719             }
720 
721           /* If quo_est was high by one, then num[i] went negative and
722              we need to correct things.  */
723           if (num[num_hi_sig] < (HOST_WIDE_INT) carry)
724             {
725               quo_est--;
726               carry = 0;                /* add divisor back in */
727               for (j = 0; j <= den_hi_sig; j++)
728                 {
729                   work = num[i + j] + den[j] + carry;
730                   carry = HIGHPART (work);
731                   num[i + j] = LOWPART (work);
732                 }
733 
734               num [num_hi_sig] += carry;
735             }
736 
737           /* Store the quotient digit.  */
738           quo[i] = quo_est;
739         }
740     }
741 


  • 商をまとめる

742   decode (quo, lquo, hquo);
743 


744  finish_up:
745   /* If result is negative, make it so.  */
746   if (quo_neg)
747     neg_double (*lquo, *hquo, lquo, hquo);
748 

  • 余りを求めるのを試す
    • このあと小数点切り上げ・切り捨てが行われる可能性がある
    • 計算式自体は最後と同じ

749   /* Compute trial remainder:  rem = num - (quo * den)  */
750   mul_double (*lquo, *hquo, lden_orig, hden_orig, lrem, hrem);
751   neg_double (*lrem, *hrem, lrem, hrem);
752   add_double (lnum_orig, hnum_orig, *lrem, *hrem, lrem, hrem);
753 


754   switch (code)
755     {
756     case TRUNC_DIV_EXPR:
757     case TRUNC_MOD_EXPR:        /* round toward zero */
758     case EXACT_DIV_EXPR:        /* for this one, it shouldn't matter */
759       return overflow;
760 

  • 負の数の切り上げ

761     case FLOOR_DIV_EXPR:
762     case FLOOR_MOD_EXPR:        /* round toward negative infinity */
763       if (quo_neg && (*lrem != 0 || *hrem != 0))   /* ratio < 0 && rem != 0 */
764         {
765           /* quo = quo - 1;  */
766           add_double (*lquo, *hquo, (HOST_WIDE_INT) -1, (HOST_WIDE_INT)  -1,
767                       lquo, hquo);
768         }
769       else
770         return overflow;
771       break;
772

  • 正の数の小数点切り上げ

773     case CEIL_DIV_EXPR:
774     case CEIL_MOD_EXPR:         /* round toward positive infinity */
775       if (!quo_neg && (*lrem != 0 || *hrem != 0))  /* ratio > 0 && rem != 0 */
776         {
777           add_double (*lquo, *hquo, (HOST_WIDE_INT) 1, (HOST_WIDE_INT) 0,
778                       lquo, hquo);
779         }
780       else
781         return overflow;
782       break;
783 

  • 四捨五入

784     case ROUND_DIV_EXPR:
785     case ROUND_MOD_EXPR:        /* round to closest integer */
786       {
787         unsigned HOST_WIDE_INT labs_rem = *lrem;
788         HOST_WIDE_INT habs_rem = *hrem;
789         unsigned HOST_WIDE_INT labs_den = lden, ltwice;
790         HOST_WIDE_INT habs_den = hden, htwice;
791 

  • 絶対値を求める

792         /* Get absolute values.  */
793         if (*hrem < 0)
794           neg_double (*lrem, *hrem, &labs_rem, &habs_rem);
795         if (hden < 0)
796           neg_double (lden, hden, &labs_den, &habs_den);
797 

  • 余りを倍にしたとき、割る数以上ならば切り上げ

798         /* If (2 * abs (lrem) >= abs (lden)) */
799         mul_double ((HOST_WIDE_INT) 2, (HOST_WIDE_INT) 0,
800                     labs_rem, habs_rem, &ltwice, &htwice);
801 
802         if (((unsigned HOST_WIDE_INT) habs_den
803              < (unsigned HOST_WIDE_INT) htwice)
804             || (((unsigned HOST_WIDE_INT) habs_den
805                  == (unsigned HOST_WIDE_INT) htwice)
806                 && (labs_den < ltwice)))
807           {
808             if (*hquo < 0)
809               /* quo = quo - 1;  */
810               add_double (*lquo, *hquo,
811                           (HOST_WIDE_INT) -1, (HOST_WIDE_INT) -1, lquo, hquo);
812             else
813               /* quo = quo + 1; */
814               add_double (*lquo, *hquo, (HOST_WIDE_INT) 1, (HOST_WIDE_INT) 0,
815                           lquo, hquo);
816           }
817         else
818           return overflow;
819       }
820       break;
821 
822     default:
823       gcc_unreachable ();
824     }
825 

  • 余りを求める

826   /* Compute true remainder:  rem = num - (quo * den)  */
827   mul_double (*lquo, *hquo, lden_orig, hden_orig, lrem, hrem);
828   neg_double (*lrem, *hrem, lrem, hrem);
829   add_double (lnum_orig, hnum_orig, *lrem, *hrem, lrem, hrem);
830   return overflow;
831 }



リンク元

Advertisement