Parent

Class Index [+]

Quicksearch

Flt::Support::Reader

Clinger algorithms to read floating point numbers from text literals with correct rounding. from his paper: “How to Read Floating Point Numbers Accurately” (William D. Clinger)

Public Class Methods

float_min_max_adj_exp(base, normalized=false) click to toggle source

Minimum & maximum adjusted exponent for numbers in base to be in the range of Floats

     # File lib/flt/support.rb, line 587
587:         def float_min_max_adj_exp(base, normalized=false)
588:           k = normalized ? base : -base
589:           unless min_max = @float_min_max_exp_values[k]
590:             max_exp = (Math.log(Float::MAX)/Math.log(base)).floor
591:             e = Float::MIN_EXP
592:             e -= Float::MANT_DIG unless normalized
593:             min_exp = (e*Math.log(Float::RADIX)/Math.log(base)).ceil
594:             @float_min_max_exp_values[k] = min_max = [min_exp, max_exp]
595:           end
596:           min_max.map{|exp| exp - 1} # adjust
597:         end
new(options={}) click to toggle source

There are two different reading approaches, selected by the :mode parameter:

  • :fixed (the destination context defines the resulting precision) input is rounded as specified by the context; if the context precision is ‘exact’, the exact input value will be represented in the destination base, which can lead to a Inexact exception (or a NaN result and an Inexact flag)

  • :free The input precision is preserved, and the destination context precision is ignored; in this case the result can be converted back to the original number (with the same precision) a rounding mode for the back conversion may be passed; otherwise any round-to-nearest is assumed. (to increase the precision of the result the input precision must be increased —adding trailing zeros)

  • :short is like :free, but the minumum number of digits that preserve the original value are generated (with :free, all significant digits are generated)

For the fixed mode there are three conversion algorithms available that can be selected with the :algorithm parameter:

  • :A Arithmetic algorithm, using correctly rounded Flt::Num arithmetic.

  • :M The Clinger Algorithm M is the slowest method, but it was the first implemented and testes and is kept as a reference for testing.

  • :R The Clinger Algorithm R, which requires an initial approximation is currently only implemented for Float and is the fastest by far.

     # File lib/flt/support.rb, line 400
400:       def initialize(options={})
401:         @exact = nil
402:         @algorithm = options[:algorithm]
403:         @mode = options[:mode] || :fixed
404:       end
ratio_float(context, u, v, k, round_mode) click to toggle source

Given exact positive integers u and v with beta**(n-1) <= u/v < beta**n and exact integer k, returns the floating point number closest to u/v * beta**n (beta is the floating-point radix)

     # File lib/flt/support.rb, line 750
750:       def self.ratio_float(context, u, v, k, round_mode)
751:         # since this handles only positive numbers and ceiling and floor
752:         # are not symmetrical, they should have been swapped before calling this.
753:         q = u.div v
754:         r = u-q*v
755:         v_r = v-r
756:         z = context.Num(1,q,k)
757:         exact = (r==0)
758:         if round_mode == :down
759:           # z = z
760:         elsif (round_mode == :up) && r>0
761:           z = context.next_plus(z)
762:         elsif r<v_r
763:           # z = z
764:         elsif r>v_r
765:           z = context.next_plus(z)
766:         else
767:           # tie
768:           if (round_mode == :half_down) || (round_mode == :half_even && ((q%2)==0)) || (round_mode == :down)
769:              # z = z
770:           else
771:             z = context.next_plus(z)
772:           end
773:         end
774:         return z, exact
775:       end

Public Instance Methods

_alg_m(context, round_mode, sign, f, e, eb, n) click to toggle source

Algorithm M to read floating point numbers from text literals with correct rounding from his paper: “How to Read Floating Point Numbers Accurately” (William D. Clinger)

     # File lib/flt/support.rb, line 709
709:       def _alg_m(context, round_mode, sign, f, e, eb, n)
710:         if e<0
711:          u,v,k = f,eb**(-e),0
712:         else
713:           u,v,k = f*(eb**e),1,0
714:         end
715:         min_e = context.etiny
716:         max_e = context.etop
717:         rp_n = context.int_radix_power(n)
718:         rp_n_1 = context.int_radix_power(n-1)
719:         r = context.radix
720:         loop do
721:            x = u.div(v) # bottleneck
722:            if (x>=rp_n_1 && x<rp_n) || k==min_e || k==max_e
723:               z, exact = Reader.ratio_float(context,u,v,k,round_mode)
724:               @exact = exact
725:               if context.respond_to?(:exception)
726:                 if k==min_e
727:                   context.exception(Num::Subnormal) if z.subnormal?
728:                   context.exception(Num::Underflow,"Input literal out of range") if z.zero? && f!=0
729:                 elsif k==max_e
730:                   if !context.exact? && z.coefficient > context.maximum_coefficient
731:                     context.exception(Num::Overflow,"Input literal out of range")
732:                   end
733:                 end
734:                 context.exception Num::Inexact if !exact
735:               end
736:               return z.copy_sign(sign)
737:            elsif x<rp_n_1
738:              u *= r
739:              k -= 1
740:            elsif x>=rp_n
741:              v *= r
742:              k += 1
743:            end
744:         end
745:       end
_alg_r(z0, context, round_mode, sign, f, e, eb, n) click to toggle source
     # File lib/flt/support.rb, line 552
552:       def _alg_r(z0, context, round_mode, sign, f, e, eb, n) # Fast for Float
553:         #raise InvalidArgument, "Reader Algorithm R only supports base 2" if context.radix != 2
554: 
555:         @z = z0
556:         @r = context.radix
557:         @rp_n_1 = context.int_radix_power(n-1)
558:         @round_mode = round_mode
559: 
560:         ret = nil
561:         loop do
562:           m, k = context.to_int_scale(@z)
563:           # TODO: replace call to compare by setting the parameters in local variables,
564:           #       then insert the body of compare here;
565:           #       then eliminate innecesary instance variables
566:           if e >= 0 && k >= 0
567:             ret = compare m, f*eb**e, m*@r**k, context
568:           elsif e >= 0 && k < 0
569:             ret = compare m, f*eb**e*@r**(-k), m, context
570:           elsif e < 0 && k >= 0
571:             ret = compare m, f, m*@r**k*eb**(-e), context
572:           else # e < 0 && k < 0
573:             ret = compare m, f*@r**(-k), m*eb**(-e), context
574:           end
575:           break if ret
576:         end
577:         ret && context.copy_sign(ret, sign) # TODO: normalize?
578:       end
_alg_r_approx(context, round_mode, sign, f, e, eb, n) click to toggle source
     # File lib/flt/support.rb, line 501
501:       def _alg_r_approx(context, round_mode, sign, f, e, eb, n)
502: 
503:         return nil if context.radix != Float::RADIX || context.exact? || context.precision > Float::MANT_DIG
504: 
505:         # Compute initial approximation; if Float uses IEEE-754 binary arithmetic, the approximation
506:         # is good enough to be adjusted in just one step.
507:         @good_approx = true
508: 
509:         ndigits = Support::AuxiliarFunctions._ndigits(f, eb)
510:         adj_exp = e + ndigits - 1
511:         min_exp, max_exp = Reader.float_min_max_adj_exp(eb)
512: 
513:         if adj_exp >= min_exp && adj_exp <= max_exp
514:           if eb==2
515:             z0 = Math.ldexp(f,e)
516:           elsif eb==10
517:             unless Flt.float_correctly_rounded?
518:               min_exp_norm, max_exp_norm = Reader.float_min_max_adj_exp(eb, true)
519:               @good_approx = false
520:               return nil if e <= min_exp_norm
521:             end
522:             z0 = Float("#{f}E#{e}")
523:           else
524:             ff = f
525:             ee = e
526:             min_exp_norm, max_exp_norm = Reader.float_min_max_adj_exp(eb, true)
527:             if e <= min_exp_norm
528:               # avoid loss of precision due to gradual underflow
529:               return nil if e <= min_exp
530:               @good_approx = false
531:               ff = Float(f)*Float(eb)**(e-min_exp_norm-1)
532:               ee = min_exp_norm + 1
533:             end
534:             # if ee < 0
535:             #   z0 = Float(ff)/Float(eb**(-ee))
536:             # else
537:             #   z0 = Float(ff)*Float(eb**ee)
538:             # end
539:             z0 = Float(ff)*Float(eb)**ee
540:           end
541: 
542:           if z0 && context.num_class != Float
543:             @good_approx = false
544:             z0 = context.Num(z0).plus(context) # context.plus(z0) ?
545:           else
546:             z0 = context.Num(z0)
547:           end
548:         end
549: 
550:       end
compare(m, x, y, context) click to toggle source
     # File lib/flt/support.rb, line 600
600:       def compare(m, x, y, context)
601:         ret = nil
602:         d = x-y
603:         d2 = 2*m*d.abs
604: 
605:         # v = f*eb**e is the number to be approximated
606:         # z = m*@r**k is the current aproximation
607:         # the error of @z is eps = abs(v-z) = 1/2 * d2 / y
608:         # we have x, y integers such that x/y = v/z
609:         # so eps < 1/2 <=> d2 < y
610:         #    d < 0 <=> x < y <=> v < z
611: 
612:         directed_rounding = [:up, :down].include?(@round_mode)
613: 
614:         if directed_rounding
615:           if @round_mode==:up ? (d <= 0) : (d < 0)
616:             # v <(=) z
617:             chk = (m == @rp_n_1) ? d2*@r : d2
618:             if (@round_mode == :up) && (chk < 2*y)
619:               # eps < 1
620:               ret = @z
621:             else
622:               @z = context.next_minus(@z)
623:             end
624:           else # @round_mode==:up ? (d > 0) : (d >= 0)
625:             # v >(=) z
626:             if (@round_mode == :down) && (d2 < 2*y)
627:               # eps < 1
628:               ret = @z
629:             else
630:               @z = context.next_plus(@z)
631:             end
632:           end
633:         else
634:           if d2 < y # eps < 1/2
635:             if (m == @rp_n_1) && (d < 0) && (y < @r*d2)
636:               # z has the minimum normalized significand, i.e. is a power of @r
637:               # and v < z
638:               # and @r*eps > 1/2
639:               # On the left of z the ulp is 1/@r than the ulp on the right; if v < z we
640:               # must require an error @r times smaller.
641:               @z = context.next_minus(@z)
642:             else
643:               # unambiguous nearest
644:               ret = @z
645:             end
646:           elsif d2 == y # eps == 1/2
647:             # round-to-nearest tie
648:             if @round_mode == :half_even
649:               if (m%2) == 0
650:                 # m is even
651:                 if (m == @rp_n_1) && (d < 0)
652:                   # z is power of @r and v < z; this wasn't really a tie because
653:                   # there are closer values on the left
654:                   @z = context.next_minus(@z)
655:                 else
656:                   # m is even => round tie to z
657:                   ret = @z
658:                 end
659:               elsif d < 0
660:                 # m is odd, v < z => round tie to prev
661:                 ret = context.next_minus(@z)
662:               elsif d > 0
663:                 # m is odd, v > z => round tie to next
664:                 ret = context.next_plus(@z)
665:               end
666:             elsif @round_mode == :half_up
667:               if d < 0
668:                 # v < z
669:                 if (m == @rp_n_1)
670:                   # this was not really a tie
671:                   @z = context.next_minus(@z)
672:                 else
673:                   ret = @z
674:                 end
675:               else # d > 0
676:                 # v >= z
677:                 ret = context.next_plus(@z)
678:               end
679:             else # @round_mode == :half_down
680:               if d < 0
681:                 # v < z
682:                 if (m == @rp_n_1)
683:                   # this was not really a tie
684:                   @z = context.next_minus(@z)
685:                 else
686:                   ret = context.next_minus(@z)
687:                 end
688:               else # d < 0
689:                 # v > z
690:                 ret = @z
691:               end
692:             end
693:           elsif d < 0 # eps > 1/2 and v < z
694:             @z = context.next_minus(@z)
695:           elsif d > 0 # eps > 1/2 and v > z
696:             @z = context.next_plus(@z)
697:           end
698:         end
699: 
700:         # Assume the initial approx is good enough (uses IEEE-754 arithmetic with round-to-nearest),
701:         # so we can avoid further iteration, except for directed rounding
702:         ret ||= @z unless directed_rounding || !@good_approx
703: 
704:         return ret
705:       end
exact?() click to toggle source
     # File lib/flt/support.rb, line 406
406:       def exact?
407:         @exact
408:       end
read(context, round_mode, sign, f, e, eb=10) click to toggle source

Given exact integers f and e, with f nonnegative, returns the floating-point number closest to f * eb**e (eb is the input radix)

If the context precision is exact an Inexact exception may occur (an NaN be returned) if an exact conversion is not possible.

round_mode: in :fixed mode it specifies how to round the result (to the context precision); it is passed separate from context for flexibility. in :free mode it specifies what rounding would be used to convert back the output to the input base eb (using the same precision that f has).

     # File lib/flt/support.rb, line 421
421:       def read(context, round_mode, sign, f, e, eb=10)
422:         @exact = true
423: 
424:         case @mode
425:         when :free, :short
426:           all_digits = (@mode == :free)
427:           # for free mode, (any) :nearest rounding is used by default
428:           Num.convert(Num[eb].Num(sign, f, e), context.num_class, :rounding=>round_mode||:nearest, :all_digits=>all_digits)
429:         when :fixed
430:           if exact_mode = context.exact?
431:             a,b = [eb, context.radix].sort
432:             m = (Math.log(b)/Math.log(a)).round
433:             if b == a**m
434:               # conmensurable bases
435:               if eb > context.radix
436:                 n = AuxiliarFunctions._ndigits(f, eb)*m
437:               else
438:                 n = (AuxiliarFunctions._ndigits(f, eb)+m-1)/m
439:               end
440:             else
441:               # inconmesurable bases; exact result may not be possible
442:               x = Num[eb].Num(sign, f, e)
443:               x = Num.convert_exact(x, context.num_class, context)
444:               @exact = !x.nan?
445:               return x
446:             end
447:           else
448:             n = context.precision
449:           end
450:           if round_mode == :nearest
451:             # :nearest is not meaningful here in :fixed mode; replace it
452:             if [:half_even, :half_up, :half_down].include?(context.rounding)
453:               round_mode = context.rounding
454:             else
455:               round_mode = :half_even
456:             end
457:           end
458:           # for fixed mode, use the context rounding by default
459:           round_mode ||= context.rounding
460:           alg = @algorithm
461:           if (context.radix == 2 && alg.nil?) || alg==:R
462:             z0 =  _alg_r_approx(context, round_mode, sign, f, e, eb, n)
463:             alg = z0 && :R
464:           end
465:           alg ||= :A
466:           case alg
467:           when :M, :R
468:             round_mode = Support.simplified_round_mode(round_mode, sign == 1)
469:             case alg
470:             when :M
471:               _alg_m(context, round_mode, sign, f, e, eb, n)
472:             when :R
473:               _alg_r(z0, context, round_mode, sign, f, e, eb, n)
474:             end
475:           else # :A
476:             # direct arithmetic conversion
477:             if round_mode == context.rounding
478:               x = Num.convert_exact(Num[eb].Num(sign, f, e), context.num_class, context)
479:               x = context.normalize(x) unless !context.respond_to?(:normalize) || context.exact?
480:               x
481:             else
482:               if context.num_class == Float
483:                 float = true
484:                 context = BinNum::FloatContext
485:               end
486:               x = context.num_class.context(context) do |local_context|
487:                 local_context.rounding = round_mode
488:                 Num.convert_exact(Num[eb].Num(sign, f, e), local_context.num_class, local_context)
489:               end
490:               if float
491:                 x = x.to_f
492:               else
493:                 x = context.normalize(x) unless context.exact?
494:               end
495:               x
496:             end
497:           end
498:         end
499:       end

Disabled; run with --debug to generate this.

[Validate]

Generated with the Darkfish Rdoc Generator 1.1.6.