in reply to Re^11: Math::BigFloat to native double?
in thread Math::BigFloat to native double?

If we limit ourselves to Math::BigFloat, I just don't see any easy way of escaping the fact that, in i2dd(), Math::BigFloat->new( sprintf "%.17f", $hi ) is simply not a good enough base 10 approximation of $hi.

For the given example, you're subtracting 3.1415926535897931000000000000000 from 3.1415926535897932384626433832795 in order to get the least significant double.
But you need to be subtracting 3.1415926535897931159979634685442 which is a far more accurate base 10 approximation of 0x1.921fb54442d18p+1. That results in a value of 0.0000000000000001224646799147353, which has a hex representation of 0x1.1a62633145c06p-53. Correctly it would be 0x1.1a62633145c07p-53 - the discrepancy being due to the rounding done in converting to base 10 values for the calculations.

In any case, the problem is that I can't see any simple way of coercing Math::BigFloat into arriving at more accurate values - as it apparently can't convert hex floats to Math::BigFloat values.

All of which brings me back to using something like mpfr.
AFAIK, my Math::MPFR ppm packages will work fine with your MSVC-built perls, but it has probably been a while since this has been properly tested.
The ppd file specifies a couple of prerequisites to install, one of which (Math::Float128) may not work correctly for you. The other is Math::Decimal64.
If you want to avoid installing those pre-requisite modules, just save a copy of my scratchpad in the cwd as Math-MPFR.ppd and run:
ppm install Math-MPFR.ppd --force
Otherwise you can run:
ppm install http://www.sisyphusion.tk/ppm/Math-Decimal64.ppd --force
ppm install http://www.sisyphusion.tk/ppm/Math-Float128.ppd --force
ppm install http://www.sisyphusion.tk/ppm/Math-MPFR.ppd --force
The prerequisites are needed only if you want to call a function that takes/returns a Math::Float128 or Math::Decimal64 object. (I don't foresee that you have any need to be calling those functions.)

I fully understand and really don't mind at all if you don't want to go down that path - and I'm sure that other alternatives exist.
But I'm pretty much out of ideas (and energy) if Math::MPFR is to be excluded. (Without it, there's just too many hoops for me).
With Math::MPFR installed, sub i2dd becomes (eg): 
sub i2dd {
  use Math::MPFR qw(:mpfr);
  my $prec = Rmpfr_get_default_prec();
  Rmpfr_set_default_prec(2098);
  my $val = Math::MPFR->new($_[0]);
  my $hi = Rmpfr_get_d($val, MPFR_RNDN);
  $val -= $hi;
  my $lo = Rmpfr_get_d($val, MPFR_RNDN);

  # Revert to original precision:
  Rmpfr_set_default_prec($prec);

  return($hi, $lo);
}

[download]
(Of course, the sub can be set up to operate more efficiently than as presented - I just wanted it to cover all angles.)

Cheers,
Rob

Replies are listed 'Best First'.
Re^13: Math::BigFloat to native double?
by syphilis (Archbishop) on Jul 15, 2015 at 14:58 UTC
    By my reckoning, this produces a correct doubledouble for 3.1415926535897932384626433832795 - and it uses only Math::BIgFloat.
    It relies on %a printf formatting, so requires perl 5.22.0. 
    use strict;
use warnings;
use Math::BigFloat;
use 5.22.0;

my %bits_table = (
 1 => '0001',
 2 => '0010',
 3 => '0011',
 4 => '0100',
 5 => '0101',
 6 => '0110',
 7 => '0111',
 8 => '1000',
 9 => '1001',
 a => '1010',
 A => '1010',
 b => '1011',
 B => '1011',
 c => '1100',
 C => '1100',
 d => '1101',
 D => '1101',
 e => '1110',
 E => '1110',
 f => '1111',
 F => '1111',
);

my %vals_table = (
0 => '1',
1 => '0.5',
2 => '0.25',
3 => '0.125',
4 => '0.0625',
5 => '0.03125',
6 => '0.015625',
7 => '0.0078125',
8 => '0.00390625',
9 => '0.001953125',
10 => '0.0009765625',
11 => '0.00048828125',
12 => '0.000244140625',
13 => '0.0001220703125',
14 => '0.00006103515625',
15 => '0.000030517578125',
16 => '0.0000152587890625',
17 => '0.00000762939453125',
18 => '0.000003814697265625',
19 => '0.0000019073486328125',
20 => '0.00000095367431640625',
21 => '0.000000476837158203125',
22 => '0.0000002384185791015625',
23 => '0.00000011920928955078125',
24 => '0.000000059604644775390625',
25 => '0.0000000298023223876953125',
26 => '0.00000001490116119384765625',
27 => '0.000000007450580596923828125',
28 => '0.0000000037252902984619140625',
29 => '0.00000000186264514923095703125',
30 => '0.000000000931322574615478515625',
31 => '0.0000000004656612873077392578125',
32 => '0.00000000023283064365386962890625',
33 => '0.000000000116415321826934814453125',
34 => '0.0000000000582076609134674072265625',
35 => '0.00000000002910383045673370361328125',
36 => '0.000000000014551915228366851806640625',
37 => '0.0000000000072759576141834259033203125',
38 => '0.00000000000363797880709171295166015625',
39 => '0.000000000001818989403545856475830078125',
40 => '0.0000000000009094947017729282379150390625',
41 => '0.00000000000045474735088646411895751953125',
42 => '0.000000000000227373675443232059478759765625',
43 => '0.0000000000001136868377216160297393798828125',
44 => '0.00000000000005684341886080801486968994140625',
45 => '0.000000000000028421709430404007434844970703125',
46 => '0.0000000000000142108547152020037174224853515625',
47 => '0.00000000000000710542735760100185871124267578125',
48 => '0.000000000000003552713678800500929355621337890625',
49 => '0.0000000000000017763568394002504646778106689453125',
50 => '0.00000000000000088817841970012523233890533447265625',
51 => '0.000000000000000444089209850062616169452667236328125',
52 => '0.0000000000000002220446049250313080847263336181640625',
53 => '0.00000000000000011102230246251565404236316680908203125',
54 => '0.000000000000000055511151231257827021181583404541015625',
55 => '0.0000000000000000277555756156289135105907917022705078125',
56 => '0.00000000000000001387778780781445675529539585113525390625',
57 => '0.000000000000000006938893903907228377647697925567626953125',
58 => '0.0000000000000000034694469519536141888238489627838134765625',
59 => '0.00000000000000000173472347597680709441192448139190673828125',
60 => '0.000000000000000000867361737988403547205962240695953369140625'
+,
61 => '0.0000000000000000004336808689942017736029811203479766845703125
+',
62 => '0.0000000000000000002168404344971008868014905601739883422851562
+5',
63 => '0.0000000000000000001084202172485504434007452800869941711425781
+25',
64 => '0.0000000000000000000542101086242752217003726400434970855712890
+625',
65 => '0.0000000000000000000271050543121376108501863200217485427856445
+3125',
66 => '0.0000000000000000000135525271560688054250931600108742713928222
+65625',
67 => '0.0000000000000000000067762635780344027125465800054371356964111
+328125',
68 => '0.0000000000000000000033881317890172013562732900027185678482055
+6640625',
69 => '0.0000000000000000000016940658945086006781366450013592839241027
+83203125',
70 => '0.0000000000000000000008470329472543003390683225006796419620513
+916015625',
71 => '0.0000000000000000000004235164736271501695341612503398209810256
+9580078125',
72 => '0.0000000000000000000002117582368135750847670806251699104905128
+4790039062',
73 => '0.0000000000000000000001058791184067875423835403125849552452564
+2395019531',
74 => '0.0000000000000000000000529395592033937711917701562924776226282
+11975097656',
75 => '0.0000000000000000000000264697796016968855958850781462388113141
+05987548828',
76 => '0.0000000000000000000000132348898008484427979425390731194056570
+52993774414',
77 => '0.0000000000000000000000066174449004242213989712695365597028285
+26496887207',
78 => '0.0000000000000000000000033087224502121106994856347682798514142
+632484436035',
79 => '0.0000000000000000000000016543612251060553497428173841399257071
+316242218018',
80 => '0.0000000000000000000000008271806125530276748714086920699628535
+6581211090088',
81 => '0.0000000000000000000000004135903062765138374357043460349814267
+8290605545044',
82 => '0.0000000000000000000000002067951531382569187178521730174907133
+9145302772522',
83 => '0.0000000000000000000000001033975765691284593589260865087453566
+9572651386261',
84 => '0.0000000000000000000000000516987882845642296794630432543726783
+47863256931305',
85 => '0.0000000000000000000000000258493941422821148397315216271863391
+73931628465652',
86 => '0.0000000000000000000000000129246970711410574198657608135931695
+86965814232826',
87 => '0.0000000000000000000000000064623485355705287099328804067965847
+934829071164131',
88 => '0.0000000000000000000000000032311742677852643549664402033982923
+967414535582066',
89 => '0.0000000000000000000000000016155871338926321774832201016991461
+983707267791033',
90 => '0.0000000000000000000000000008077935669463160887416100508495730
+9918536338955164',
91 => '0.0000000000000000000000000004038967834731580443708050254247865
+4959268169477582',
92 => '0.0000000000000000000000000002019483917365790221854025127123932
+7479634084738791',
93 => '0.0000000000000000000000000001009741958682895110927012563561966
+3739817042369395',
94 => '0.0000000000000000000000000000504870979341447555463506281780983
+18699085211846977',
95 => '0.0000000000000000000000000000252435489670723777731753140890491
+59349542605923489',
96 => '0.0000000000000000000000000000126217744835361888865876570445245
+79674771302961744',
97 => '0.0000000000000000000000000000063108872417680944432938285222622
+898373856514808722',
98 => '0.0000000000000000000000000000031554436208840472216469142611311
+449186928257404361',
99 => '0.0000000000000000000000000000015777218104420236108234571305655
+72459346412870218',
100 => '0.000000000000000000000000000000788860905221011805411728565282
+78622967320643510902',
101 => '0.000000000000000000000000000000394430452610505902705864282641
+39311483660321755451',
102 => '0.000000000000000000000000000000197215226305252951352932141320
+69655741830160877726',
103 => '0.000000000000000000000000000000098607613152626475676466070660
+348278709150804388628',
104 => '0.000000000000000000000000000000049303806576313237838233035330
+174139354575402194314',
105 => '0.000000000000000000000000000000024651903288156618919116517665
+087069677287701097157',
106 => '0.000000000000000000000000000000012325951644078309459558258832
+543534838643850548578',
107 => '0.000000000000000000000000000000006162975822039154729779129416
+2717674193219252742892',
108 => '0.000000000000000000000000000000003081487911019577364889564708
+1358837096609626371446',
109 => '0.000000000000000000000000000000001540743955509788682444782354
+0679418548304813185723',
110 => '0.000000000000000000000000000000000770371977754894341222391177
+03397092741524065928616',
111 => '0.000000000000000000000000000000000385185988877447170611195588
+51698546370762032964308',
112 => '0.000000000000000000000000000000000192592994438723585305597794
+25849273185381016482154',
113 => '0.000000000000000000000000000000000096296497219361792652798897
+129246365926905082410769',
114 => '0.000000000000000000000000000000000048148248609680896326399448
+564623182963452541205385',
115 => '0.000000000000000000000000000000000024074124304840448163199724
+282311591481726270602692',
116 => '0.000000000000000000000000000000000012037062152420224081599862
+141155795740863135301346',
117 => '0.000000000000000000000000000000000006018531076210112040799931
+0705778978704315676506731',
118 => '0.000000000000000000000000000000000003009265538105056020399965
+5352889489352157838253365',
119 => '0.000000000000000000000000000000000001504632769052528010199982
+7676444744676078919126683',
120 => '0.000000000000000000000000000000000000752316384526264005099991
+38382223723380394595633414',
121 => '0.000000000000000000000000000000000000376158192263132002549995
+69191111861690197297816707',
122 => '0.000000000000000000000000000000000000188079096131566001274997
+84595555930845098648908353',
123 => '0.000000000000000000000000000000000000094039548065783000637498
+922977779654225493244541767',
124 => '0.000000000000000000000000000000000000047019774032891500318749
+461488889827112746622270884',
125 => '0.000000000000000000000000000000000000023509887016445750159374
+730744444913556373311135442',
126 => '0.000000000000000000000000000000000000011754943508222875079687
+365372222456778186655567721',
127 => '0.000000000000000000000000000000000000005877471754111437539843
+6826861112283890933277838604',
128 => '0.000000000000000000000000000000000000002938735877055718769921
+8413430556141945466638919302',
129 => '0.000000000000000000000000000000000000001469367938527859384960
+9206715278070972733319459651',
130 => '0.000000000000000000000000000000000000000734683969263929692480
+46033576390354863666597298255',
131 => '0.000000000000000000000000000000000000000367341984631964846240
+23016788195177431833298649128',
132 => '0.000000000000000000000000000000000000000183670992315982423120
+11508394097588715916649324564',
133 => '0.000000000000000000000000000000000000000091835496157991211560
+057541970487943579583246622819',
134 => '0.000000000000000000000000000000000000000045917748078995605780
+02877098524397178979162331141',
135 => '0.000000000000000000000000000000000000000022958874039497802890
+014385492621985894895811655705',
136 => '0.000000000000000000000000000000000000000011479437019748901445
+007192746310992947447905827852',
137 => '0.000000000000000000000000000000000000000005739718509874450722
+5035963731554964737239529139262',
138 => '0.000000000000000000000000000000000000000002869859254937225361
+2517981865777482368619764569631',
139 => '0.000000000000000000000000000000000000000001434929627468612680
+6258990932888741184309882284816',
140 => '0.000000000000000000000000000000000000000000717464813734306340
+31294954664443705921549411424078',
141 => '0.000000000000000000000000000000000000000000358732406867153170
+15647477332221852960774705712039',
142 => '0.000000000000000000000000000000000000000000179366203433576585
+07823738666110926480387352856019',
143 => '0.000000000000000000000000000000000000000000089683101716788292
+539118693330554632401936764280097',
144 => '0.000000000000000000000000000000000000000000044841550858394146
+269559346665277316200968382140049',
145 => '0.000000000000000000000000000000000000000000022420775429197073
+134779673332638658100484191070024',
146 => '0.000000000000000000000000000000000000000000011210387714598536
+567389836666319329050242095535012',
147 => '0.000000000000000000000000000000000000000000005605193857299268
+2836949183331596645251210477675061',
148 => '0.000000000000000000000000000000000000000000002802596928649634
+141847459166579832262560523883753',
149 => '0.000000000000000000000000000000000000000000001401298464324817
+0709237295832899161312802619418765',
150 => '0.000000000000000000000000000000000000000000000700649232162408
+53546186479164495806564013097093826',
151 => '0.000000000000000000000000000000000000000000000350324616081204
+26773093239582247903282006548546913',
152 => '0.000000000000000000000000000000000000000000000175162308040602
+13386546619791123951641003274273456',
153 => '0.000000000000000000000000000000000000000000000087581154020301
+066932733098955619758205016371367282',
154 => '0.000000000000000000000000000000000000000000000043790577010150
+533466366549477809879102508185683641',
155 => '0.000000000000000000000000000000000000000000000021895288505075
+266733183274738904939551254092841821',
156 => '0.000000000000000000000000000000000000000000000010947644252537
+63336659163736945246977562704642091',
157 => '0.000000000000000000000000000000000000000000000005473822126268
+8166832958186847262348878135232104551',
158 => '0.000000000000000000000000000000000000000000000002736911063134
+4083416479093423631174439067616052276',
159 => '0.000000000000000000000000000000000000000000000001368455531567
+2041708239546711815587219533808026138',
);

my $pi = '3.1415926535897932384626433832795';
my ($msd, $lsd) = i2dd($pi, \%bits_table, \%vals_table);
print "$msd $lsd\n";
printf "%a %a\n", $msd, $lsd ;


sub i2dd {
  my $str = $_[0];
  my $sign = $str =~ s/^\-// ? -1 : 1;
  my $val = Math::BigFloat->new($str);
  my $msd = 0 + $val->bstr;
  my ($bin, $exp) = d2bin($msd, $_[1]);
  my $to_subtract = bin2val($bin, $_[2]) * (2 ** $exp);
  $val -= $to_subtract;
  my $lsd = 0 + $val->bstr;
  return ($msd * $sign, $lsd * $sign);
}


sub d2bin {
  my ($mant, $exp) = split /p/, sprintf("%a", $_[0]);
  my %bits = %{$_[1]};
  $mant =~ s/-?0x//i;
  my $bin = $mant =~ /^1/ ? 1 : 0;
  $mant =~ s/.+\.//;
  for(1 .. length($mant)) {
    $bin .= $bits{substr($mant, $_ - 1, 1)};
  }
  return ($bin, $exp);
}

sub bin2val {
  my $bin = $_[0];
  my %vals = %{$_[1]};
  my $val = Math::BigFloat->new(0);

  for(1 .. length($bin)) {
    $val += Math::BigFloat->new($vals{$_ - 1}) if substr($bin, $_ - 1,
+ 1);  
  }
  return $val;
}

__END__

Outputs:
3.14159265358979 1.22464679914735e-016
0x1.921fb54442d18p+1 0x1.1a62633145c07p-53

    [download]
    The values in %vals_tab have been calculated to 160 bits of precision. 160 is an arbitrarily chosen value, == 1 + (3 * 53).
    I make no claims about this approach. (I'll scrutinise it and test it over the next week or so.)

    Cheers,
    Rob
[reply]
[d/l]

      Thanks syphilis. That's beyond the call.

      I'm still digesting Re^12: Math::BigFloat to native double?; but I took "But I'm pretty much out of ideas (and energy) if Math::MPFR is to be excluded. (Without it, there's just too many hoops for me)." at face value, so this is most unexpected :) Thank you.

      "(I'll scrutinise it and test it over the next week or so.)" But, please don't expend more of your energies on this. It's just a background project for me; something I do whilst I'm waiting for other stuff to run; or when I hit a mental block in my current work and need a diversion to allow my brain to unencumber itself of the details.

      Given where I started with the simple question in the OP, this thread's has gone way beyond anything I could have expected anyone else to give a damn about :)

      With the rise and rise of 'Social' network sites: 'Computers are making people easier to use everyday'
      Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
      "Science is about questioning the status quo. Questioning authority".
      In the absence of evidence, opinion is indistinguishable from prejudice.
      I'm with torvalds on this Agile (and TDD) debunked I told'em LLVM was the way to go. But did they listen!
[reply]
        Thing is that, with re-charged batteries, I start thinking "Hang on, I *can* improve on that". (And there's still room for improvement.)

        Anyway - good luck, and enjoy the exercise. I find the doubledouble to be rather intriguing (not to mention a little exasperating), and I purchased the (second hand) PPC box for the sole purpose of looking more readily at it.
        You should be able to pick one up over there for next-to-nothing if you ever get the urge. I had one offered to me for free - it was a genuine offer, too, but would've required a trip to New York to pick it up.
        Being big-endian, the PPC box also provides a handy reference to that architecture, too.

        Cheers,
        Rob
[reply]

In Section Seekers of Perl Wisdom