Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.forth > #57543 > unrolled thread

real to rational number

Back to article view | Back to comp.lang.forth

Contents

  real to rational number krishna.myneni@ccreweb.org - 2017-04-18 04:49 -0700
    Re: real to rational number lehs <skydda.barnen@gmail.com> - 2017-04-18 06:48 -0700
      Re: real to rational number krishna.myneni@ccreweb.org - 2017-04-18 19:09 -0700
        Re: real to rational number lehs <skydda.barnen@gmail.com> - 2017-04-18 22:53 -0700
          Re: real to rational number krishna.myneni@ccreweb.org - 2017-04-19 04:43 -0700
            Re: real to rational number rickman <gnuarm@gmail.com> - 2017-04-19 15:41 -0400
              Re: real to rational number Paul Rubin <no.email@nospam.invalid> - 2017-04-19 13:03 -0700
                Re: real to rational number "Rudy Velthuis" <newsgroups@rvelthuis.de> - 2017-04-20 17:26 +0000
                Re: real to rational number hughaguilar96@gmail.com - 2017-04-20 14:29 -0700
        Re: real to rational number "Rudy Velthuis" <newsgroups@rvelthuis.de> - 2017-04-19 23:44 +0000
          Re: real to rational number Paul Rubin <no.email@nospam.invalid> - 2017-04-19 19:30 -0700
            Re: real to rational number "Rudy Velthuis" <newsgroups@rvelthuis.de> - 2017-04-20 17:08 +0000
              Re: real to rational number Paul Rubin <no.email@nospam.invalid> - 2017-04-20 11:46 -0700
                Re: real to rational number "Rudy Velthuis" <newsgroups@rvelthuis.de> - 2017-04-21 09:47 +0000
              Re: real to rational number lehs <skydda.barnen@gmail.com> - 2017-04-20 20:49 -0700
      Re: real to rational number rickman <gnuarm@gmail.com> - 2017-04-19 15:36 -0400
        Re: real to rational number krishna.myneni@ccreweb.org - 2017-04-20 18:44 -0700
    Re: real to rational number hughaguilar96@gmail.com - 2017-04-18 07:06 -0700
      Re: real to rational number krishna.myneni@ccreweb.org - 2017-04-18 19:17 -0700
        Re: real to rational number Paul Rubin <no.email@nospam.invalid> - 2017-04-18 19:45 -0700
          Re: real to rational number krishna.myneni@ccreweb.org - 2017-04-18 19:57 -0700
          Re: real to rational number anton@mips.complang.tuwien.ac.at (Anton Ertl) - 2017-04-19 05:35 +0000
        Re: real to rational number hughaguilar96@gmail.com - 2017-04-18 23:42 -0700
          Re: real to rational number krishna.myneni@ccreweb.org - 2017-04-19 04:53 -0700
            Re: real to rational number hughaguilar96@gmail.com - 2017-04-19 11:32 -0700
    Re: real to rational number minforth@arcor.de - 2017-04-18 08:23 -0700
      Re: real to rational number krishna.myneni@ccreweb.org - 2017-04-18 19:24 -0700
        Re: real to rational number Paul Rubin <no.email@nospam.invalid> - 2017-04-18 19:48 -0700
          Re: real to rational number krishna.myneni@ccreweb.org - 2017-04-18 20:01 -0700
    Re: real to rational number Paul Rubin <no.email@nospam.invalid> - 2017-04-18 19:39 -0700
      Re: real to rational number krishna.myneni@ccreweb.org - 2017-04-18 19:56 -0700
        Re: real to rational number Paul Rubin <no.email@nospam.invalid> - 2017-04-19 00:41 -0700
          Re: real to rational number krishna.myneni@ccreweb.org - 2017-04-19 04:48 -0700
          Re: real to rational number hughaguilar96@gmail.com - 2017-04-19 11:30 -0700
        Re: real to rational number albert@cherry.spenarnc.xs4all.nl (Albert van der Horst) - 2017-04-20 20:15 +0200
          Re: real to rational number krishna.myneni@ccreweb.org - 2017-04-20 15:23 -0700
            Re: real to rational number Paul Rubin <no.email@nospam.invalid> - 2017-04-20 18:52 -0700
    Re: real to rational number Pablo Hugo Reda <pabloreda@gmail.com> - 2017-04-18 19:58 -0700
      Re: real to rational number krishna.myneni@ccreweb.org - 2017-04-18 20:05 -0700
        Re: real to rational number Pablo Hugo Reda <pabloreda@gmail.com> - 2017-04-19 05:30 -0700
          Re: real to rational number Pablo Hugo Reda <pabloreda@gmail.com> - 2017-04-19 05:35 -0700
            Re: real to rational number Pablo Hugo Reda <pabloreda@gmail.com> - 2017-04-19 05:41 -0700
              Re: real to rational number krishna.myneni@ccreweb.org - 2017-04-19 19:39 -0700
    Re: real to rational number "Rudy Velthuis" <newsgroups@rvelthuis.de> - 2017-04-19 23:16 +0000
      Re: real to rational number krishna.myneni@ccreweb.org - 2017-04-19 19:42 -0700
        Re: real to rational number "Rudy Velthuis" <newsgroups@rvelthuis.de> - 2017-04-20 17:15 +0000
          Re: real to rational number krishna.myneni@ccreweb.org - 2017-04-20 15:11 -0700
            Re: real to rational number hughaguilar96@gmail.com - 2017-04-20 15:49 -0700
              Re: real to rational number Alex <alex@rivadpm.com> - 2017-04-21 00:20 +0100
                Re: real to rational number krishna.myneni@ccreweb.org - 2017-04-20 18:49 -0700
                  Re: real to rational number The Beez <the.beez.speaks@gmail.com> - 2017-04-21 03:18 -0700
                    Re: real to rational number krishna.myneni@ccreweb.org - 2017-04-21 18:39 -0700
                      Re: real to rational number krishna.myneni@ccreweb.org - 2017-04-21 18:46 -0700
                      Re: real to rational number Paul Rubin <no.email@nospam.invalid> - 2017-04-21 21:55 -0700
                        Re: real to rational number krishna.myneni@ccreweb.org - 2017-04-22 08:33 -0700
                          Re: real to rational number Paul Rubin <no.email@nospam.invalid> - 2017-04-22 13:35 -0700
                      Re: real to rational number The Beez <the.beez.speaks@gmail.com> - 2017-04-24 00:58 -0700
                        Re: real to rational number minforth@arcor.de - 2017-04-24 01:08 -0700
                    Re: real to rational number minforth@arcor.de - 2017-04-23 02:28 -0700
                  Re: real to rational number hughaguilar96@gmail.com - 2017-04-21 18:43 -0700
                    rationalising  reality,   was Re: real to rational number Jan Coombs <jenfhaomndgfwutc@murmic.plus.com> - 2017-04-22 10:28 +0100
                      Re: rationalising  reality,   was Re: real to rational number Mark Wills <markwills1970@gmail.com> - 2017-04-22 02:33 -0700
              Re: real to rational number krishna.myneni@ccreweb.org - 2017-04-20 16:22 -0700
                Re: real to rational number hughaguilar96@gmail.com - 2017-04-20 22:02 -0700
                  Re: real to rational number Julian Fondren <julian.fondren@gmail.com> - 2017-04-21 00:12 -0700
                  Re: real to rational number "Rudy Velthuis" <newsgroups@rvelthuis.de> - 2017-04-21 09:55 +0000
                    Re: real to rational number hughaguilar96@gmail.com - 2017-04-25 10:52 -0700
                      Re: real to rational number "Rudy Velthuis" <newsgroups@rvelthuis.de> - 2017-04-29 03:16 +0000
                Re: real to rational number "Rudy Velthuis" <newsgroups@rvelthuis.de> - 2017-04-21 09:52 +0000
                Re: real to rational number lehs <skydda.barnen@gmail.com> - 2017-04-21 05:35 -0700
                  Re: real to rational number hughaguilar96@gmail.com - 2017-04-21 19:05 -0700
                    Re: real to rational number The Beez <the.beez.speaks@gmail.com> - 2017-04-24 01:14 -0700
                      Re: real to rational number Cecil Bayona <cbayona@cbayona.com> - 2017-04-24 11:31 -0500
                        Re: real to rational number The Beez <the.beez.speaks@gmail.com> - 2017-04-25 00:52 -0700
                      Re: real to rational number hughaguilar96@gmail.com - 2017-04-24 11:43 -0700
                        Re: real to rational number The Beez <the.beez.speaks@gmail.com> - 2017-04-25 00:35 -0700
                        Re: real to rational number The Beez <the.beez.speaks@gmail.com> - 2017-04-25 01:09 -0700
                          Re: real to rational number minforth@arcor.de - 2017-04-26 01:55 -0700

Page 2 of 4 — ← Prev page 1 [2] 3 4  Next page →

#57574

On Tuesday, April 18, 2017 at 9:45:41 PM UTC-5, Paul Rubin wrote:
> krishna.myneni@ccreweb.org writes:
> > r = 7.61904761904762e-02
> > Find the fraction a/b, so that it agrees to 15 significant digits.
> 
> Wait, that number is 8/105.  That was obvious on looking at 1/r.

Yes, indeed. Please encode your algorithm ;-)

KM

[toc] | [prev] | [next] | [standalone]

#57579

Paul Rubin <no.email@nospam.invalid> writes:
>krishna.myneni@ccreweb.org writes:
>> r = 7.61904761904762e-02
>> Find the fraction a/b, so that it agrees to 15 significant digits.
>
>Wait, that number is 8/105.  That was obvious on looking at 1/r.

8/105=.076|190476...

where the part after | repeats itself.  An exact fraction for r above is

761904761904762/100

which can be simplified by dividing by the GCD to

380952380952381/50

If you allow inexact fractions (I don't see that Rosettacode does),
things get tricky.

- anton
-- 
M. Anton Ertl  http://www.complang.tuwien.ac.at/anton/home.html
comp.lang.forth FAQs: http://www.complang.tuwien.ac.at/forth/faq/toc.html
     New standard: http://www.forth200x.org/forth200x.html
   EuroForth 2016: http://www.euroforth.org/ef16/

[toc] | [prev] | [next] | [standalone]

#57583

On Tuesday, April 18, 2017 at 7:17:07 PM UTC-7, krishna...@ccreweb.org wrote:
> On Tuesday, April 18, 2017 at 9:06:39 AM UTC-5, hughag...@gmail.com wrote:
> > I have rational numbers in the novice-package. This goes way back to 2010.
> > 
> > They don't work very well because the numerator and denominator are 32-bit, which is pretty small. They'll work a lot better in Straight Forth that requires a 64-bit cell and provides a double-stack that is separate from the single-stack.
> > 
> > I don't know what you mean by a "brute force search." I wrote code that determines the GCD --- isn't that the way that it is always done?
> 
> Take, for example, the following real number:
> 
> r = 7.61904761904762e-02
> 
> Find the fraction a/b, so that it agrees to 15 significant digits. It is a simple fraction. Can a double length GCD can do this in a reasonable time?
> 
> Krishna

When I said that I have code for GCD, I meant that this is used for adding ratios together --- they have to have a common denominator so you can add the numerators.

For converting a real number into a ratio, you need to use continued fractions. I have this in the novice-package, but I didn't write it --- it was written by Nathaniel Grossman and published in "Forth Dimensions" (Sep/Oct 1984) --- all I did was upgrade his program to work under 32-bit ANS-Forth rather than 16-bit Forth-83.

CF is very slow because ANS-Forth lacks a double-precision division --- this had to be written in high-level Forth --- this is UD/MOD which is in the support code.

Here is an example:

00761904761904762.  ok-2 
10000000000000000.  ok-4 
cf 
partial-quot             numerator           denominator
           0                     0                     1
          13                     1                    13
           8                     8                   105
9523809523809        76190476190473       999999999999958
           2       152380952380954      2000000000000021
           2       380952380952381      5000000000000000 ok

Here is another example:

3141592653589793238.  ok-2 
1000000000000000000.  ok-4 
cf 
partial-quot             numerator           denominator
           3                     3                     1
           7                    22                     7
          15                   333                   106
           1                   355                   113
         292                103993                 33102
           1                104348                 33215
           1                208341                 66317
           1                312689                 99532
           2                833719                265381
           1               1146408                364913
           3               4272943               1360120
           1               5419351               1725033
          14              80143857              25510582
           2             165707065              52746197
           1             245850922              78256779
           1             411557987             131002976
           2            1068966896             340262731
           2            2549491779             811528438
           1            3618458675            1151791169
          12           45970995879           14633022466
           1           49589454554           15784813635
           3          194739359541           61987463371
           1          244328814095           77772277006
           6         1660712244111          528621125407
           1         1905041058206          606393402413
           4         9280876476935         2954194735059
           1        11185917535141         3560588137472
          18       210627392109473        67044781209555
           1       221813309644614        70605369347027
           3       876067321043315       278860889250636
           5      4602149914861189      1464909815600207
           2     10080367150765693      3208680520451050
           2     24762884216392575      7882270856502307
           1     34843251367158268     11090951376953357
           1     59606135583550843     18973222233455664
           1     94449386950709111     30064173610409021
          16   1570796326794896619    500000000000000000 ok

I still don't know what you mean by a "brute force search."

[toc] | [prev] | [next] | [standalone]

#57594

On Wednesday, April 19, 2017 at 1:42:55 AM UTC-5, hughag...@gmail.com wrote:
> On Tuesday, April 18, 2017 at 7:17:07 PM UTC-7, krishna...@ccreweb.org wrote:
> > On Tuesday, April 18, 2017 at 9:06:39 AM UTC-5, hughag...@gmail.com wrote:
> > > I have rational numbers in the novice-package. This goes way back to 2010.
> > > 
> > > They don't work very well because the numerator and denominator are 32-bit, which is pretty small. They'll work a lot better in Straight Forth that requires a 64-bit cell and provides a double-stack that is separate from the single-stack.
> > > 
> > > I don't know what you mean by a "brute force search." I wrote code that determines the GCD --- isn't that the way that it is always done?
> > 
> > Take, for example, the following real number:
> > 
> > r = 7.61904761904762e-02
> > 
> > Find the fraction a/b, so that it agrees to 15 significant digits. It is a simple fraction. Can a double length GCD can do this in a reasonable time?
> > 
> > Krishna
> 
> When I said that I have code for GCD, I meant that this is used for adding ratios together --- they have to have a common denominator so you can add the numerators.
> 
> For converting a real number into a ratio, you need to use continued fractions. I have this in the novice-package, but I didn't write it --- it was written by Nathaniel Grossman and published in "Forth Dimensions" (Sep/Oct 1984) --- all I did was upgrade his program to work under 32-bit ANS-Forth rather than 16-bit Forth-83.
> 
> CF is very slow because ANS-Forth lacks a double-precision division --- this had to be written in high-level Forth --- this is UD/MOD which is in the support code.
> 
> Here is an example:
> 
> 00761904761904762.  ok-2 
> 10000000000000000.  ok-4 
> cf 
> partial-quot             numerator           denominator
>            0                     0                     1
>           13                     1                    13
>            8                     8                   105
> ...

Very good. The continued fraction should truncate when we have achieve the same number of significant digits specified for r, in this case, 15. With regard to your question about what constitutes a brute-force search, please see my reply to lehs.

Krishna

[toc] | [prev] | [next] | [standalone]

#57602

On Wednesday, April 19, 2017 at 4:53:13 AM UTC-7, krishna...@ccreweb.org wrote:
> On Wednesday, April 19, 2017 at 1:42:55 AM UTC-5, hughag...@gmail.com wrote:
> > On Tuesday, April 18, 2017 at 7:17:07 PM UTC-7, krishna...@ccreweb.org wrote:
> > > On Tuesday, April 18, 2017 at 9:06:39 AM UTC-5, hughag...@gmail.com wrote:
> > > > I have rational numbers in the novice-package. This goes way back to 2010.
> > > > 
> > > > They don't work very well because the numerator and denominator are 32-bit, which is pretty small. They'll work a lot better in Straight Forth that requires a 64-bit cell and provides a double-stack that is separate from the single-stack.
> > > > 
> > > > I don't know what you mean by a "brute force search." I wrote code that determines the GCD --- isn't that the way that it is always done?
> > > 
> > > Take, for example, the following real number:
> > > 
> > > r = 7.61904761904762e-02
> > > 
> > > Find the fraction a/b, so that it agrees to 15 significant digits. It is a simple fraction. Can a double length GCD can do this in a reasonable time?
> > > 
> > > Krishna
> > 
> > When I said that I have code for GCD, I meant that this is used for adding ratios together --- they have to have a common denominator so you can add the numerators.
> > 
> > For converting a real number into a ratio, you need to use continued fractions. I have this in the novice-package, but I didn't write it --- it was written by Nathaniel Grossman and published in "Forth Dimensions" (Sep/Oct 1984) --- all I did was upgrade his program to work under 32-bit ANS-Forth rather than 16-bit Forth-83.
> > 
> > CF is very slow because ANS-Forth lacks a double-precision division --- this had to be written in high-level Forth --- this is UD/MOD which is in the support code.
> > 
> > Here is an example:
> > 
> > 00761904761904762.  ok-2 
> > 10000000000000000.  ok-4 
> > cf 
> > partial-quot             numerator           denominator
> >            0                     0                     1
> >           13                     1                    13
> >            8                     8                   105
> > ...
> 
> Very good. The continued fraction should truncate when we have achieve the same number of significant digits specified for r, in this case, 15. With regard to your question about what constitutes a brute-force search, please see my reply to lehs.
> 
> Krishna

You are right. I should make it stop after it has achieved maximum precision.

[toc] | [prev] | [next] | [standalone]

#57550

Am Dienstag, 18. April 2017 13:49:31 UTC+2 schrieb krishna...@ccreweb.org:
> Rosetta Code has a number of entries for code to converting a (computer representation) of a real number into a rational number, or to the nearest expressible fraction, a/b, where both a and b are integers. There is even a Forth entry!
> 
> https://rosettacode.org/wiki/Convert_decimal_number_to_rational
> 
> However, it seems like many of the algorithms expressed there are brute force searches, and don't make use of the fundamental theorem of arithmetic to find the suitable integers a and b (though the code which uses library functions may do so).
> 
> This particular problem has been in the back of my mind for some time, and it is one that comes up from time to time. For example, in the quantum physics of angular momentum, the Wigner 3-j symbols are real numbers between zero and one, expressible as +/-sqrt(a/b), and if one computes the quantities as real numbers, it would be nice to find (a/b) in an efficient manner. In tables of these symbols, the rational number a/b is expressed as a series of powers of prime numbers and a sign term, e.g. -1 [1 -3 0 -1] represents 
> 
> r = -sqrt(2/(3^3*7)).
> 
> When the finite representation of the real value is given, say to double precision, a smart algorithm to find a and b, or the prime powers representation would be handy. The real to rational code is certainly useful, but I think the algorithms should be able to take advantage of the fact that |r| <= 1. It's an interesting inversion problem.
> 
> Incidentally, here is Forth code to compute Wigner 3j symbols -- it is somewhat limited in the range of usable angular momentum quantum numbers, though.
> 
> ftp://ccreweb.org/software/kforth/examples/fsl/extras/cg.4th
> 
> and accompanying test code,
> 
> ftp://ccreweb.org/software/kforth/examples/fsl/extras/cg-test.4th
> 
> Krishna

Rational fraction approximation for math students in Forth:

https://groups.google.com/forum/#!searchin/comp.lang.forth/rational$20fraction$20representation$20integer$20ratio%7Csort:relevance/comp.lang.forth/CS-ZKjffwgk/wuYthOQe6SUJ

[toc] | [prev] | [next] | [standalone]

#57567

On Tuesday, April 18, 2017 at 10:23:27 AM UTC-5, minf...@arcor.de wrote:
...
> Rational fraction approximation for math students in Forth:
> 
> https://groups.google.com/forum/#!searchin/comp.lang.forth/rational$20fraction$20representation$20integer$20ratio%7Csort:relevance/comp.lang.forth/CS-ZKjffwgk/wuYthOQe6SUJ

I'll check it out. Thanks.

Krishna

[toc] | [prev] | [next] | [standalone]

#57572

krishna.myneni@ccreweb.org writes:
> I'll check it out. Thanks.

Do you know about the inverse symbolic calculator?

  https://isc.carma.newcastle.edu.au/

You might like to check that out too.

[toc] | [prev] | [next] | [standalone]

#57576

On Tuesday, April 18, 2017 at 9:49:00 PM UTC-5, Paul Rubin wrote:
> krishna.myneni@ccreweb.org writes:
> > I'll check it out. Thanks.
> 
> Do you know about the inverse symbolic calculator?
> 
>   https://isc.carma.newcastle.edu.au/
> 
> You might like to check that out too.

Thanks. Wolfram Alpha will also do this, e.g.

Rationalize(7.61904761904762e-02)

will return 8/105 .

[toc] | [prev] | [next] | [standalone]

#57570

krishna.myneni@ccreweb.org writes:
> Rosetta Code has a number of entries for code to converting a
> (computer representation) of a real number into a rational number, or
> to the nearest expressible fraction, a/b, where both a and b are
> integers. There is even a Forth entry!

The Forth entry looks pretty naive unfortunately.  I'm not sure I
totally understood the problem but I thought the sequence of
lowest-terms rational approximations to a real was given by the
continued fraction sequence.

> r = -sqrt(2/(3^3*7)).
> When the finite representation of the real value is given, say to
> double precision, a smart algorithm to find a and b, or the prime
> powers representation would be handy. 

If for some reason you didn't want to do the calculation directly, like
you're on an MCU or something, maybe you could keep a table of the
logarithms of a list of primes.  Then you can find -log(r) by adding
suitable multiples of the needed primes, and dividing by 2.

[toc] | [prev] | [next] | [standalone]

#57573

On Tuesday, April 18, 2017 at 9:40:04 PM UTC-5, Paul Rubin wrote:
> krishna.myneni@ccreweb.org writes:
> > Rosetta Code has a number of entries for code to converting a
> > (computer representation) of a real number into a rational number, or
> > to the nearest expressible fraction, a/b, where both a and b are
> > integers. There is even a Forth entry!
> 
> The Forth entry looks pretty naive unfortunately. ...

Most of the entries in other languages look naive also, in the sense that they appear to be brute force searches.

> I'm not sure I
> totally understood the problem but I thought the sequence of
> lowest-terms rational approximations to a real was given by the
> continued fraction sequence.
> 
> > r = -sqrt(2/(3^3*7)).
> > When the finite representation of the real value is given, say to
> > double precision, a smart algorithm to find a and b, or the prime
> > powers representation would be handy. 
> 
> If for some reason you didn't want to do the calculation directly, like
> you're on an MCU or something, maybe you could keep a table of the
> logarithms of a list of primes.  Then you can find -log(r) by adding
> suitable multiples of the needed primes, and dividing by 2.

Just so the problem is clear, I have a double precision value for r, which means I know it to about 16 significant digits. I know r, with infinite precision, is a rational number, but my machine representation is a rounded version to the limited precision. Now, I want to find "a" and "b" where r = a/b. Dividing the floating point representations of a and b, ra and rb, should give me a number equal to r in the machine representation.

Krishna
  

[toc] | [prev] | [next] | [standalone]

#57584

krishna.myneni@ccreweb.org writes:
> Just so the problem is clear, I have a double precision value for r,
> which means I know it to about 16 significant digits. I know r, with
> infinite precision, is a rational number, but my machine
> representation is a rounded version to the limited precision. Now, I
> want to find "a" and "b" where r = a/b. Dividing the floating point
> representations of a and b, ra and rb, should give me a number equal
> to r in the machine representation.

For the r you gave, it's 8/105 in 64-bit gforth:

 7.61904761904762e-02 fvariable r  
 8e 105e f/ fvariable q
 q r f- f. 

=>   0.  ok

If you want to find the ratio by continued fractions, here's Haskell
code that I just wrote:

    import Data.Ratio
    import Control.Monad (forM_)

    r = 7.61904761904762e-02

    -- turn real number into continued fraction sequence
    contFrac :: RealFrac a => a -> [Integer]
    contFrac x = a : contFrac (1/b) where
      (a,b) = properFraction x

    -- turn truncated continued fraction into a rational number
    unFrac :: Integral a => [a] -> Ratio Integer
    unFrac [c] = fromIntegral c
    unFrac (c:cs) = fromIntegral c + 1/unFrac cs

    -- make rational from first n terms of cont frac of a real
    approx :: RealFrac a => Int -> a -> Ratio Integer
    approx n x = unFrac . take n . contFrac $ x

    -- print first 4 rational approximations of r
    main = forM_ [1..4] $ \i -> do
              print (approx i r)

The output for the first 4 terms is:

    *Main> main
    0 % 1
    1 % 13
    8 % 105
    70368744177665 % 923589767331853

You can see the 4th approximant is basically noise.  The 5th is much
worse noise, several dozen digits.

[toc] | [prev] | [next] | [standalone]

#57593

On Wednesday, April 19, 2017 at 2:41:10 AM UTC-5, Paul Rubin wrote:
> krishna.myneni@ccreweb.org writes:
> > Just so the problem is clear, I have a double precision value for r,
> > which means I know it to about 16 significant digits. I know r, with
> > infinite precision, is a rational number, but my machine
> > representation is a rounded version to the limited precision. Now, I
> > want to find "a" and "b" where r = a/b. Dividing the floating point
> > representations of a and b, ra and rb, should give me a number equal
> > to r in the machine representation.
> 
> For the r you gave, it's 8/105 in 64-bit gforth:
> 
>  7.61904761904762e-02 fvariable r  
>  8e 105e f/ fvariable q
>  q r f- f. 
> 
> =>   0.  ok
> 
> If you want to find the ratio by continued fractions, here's Haskell
> code that I just wrote:
> 
>     import Data.Ratio
>     import Control.Monad (forM_)
> 
>     r = 7.61904761904762e-02
> 
>     -- turn real number into continued fraction sequence
>     contFrac :: RealFrac a => a -> [Integer]
>     contFrac x = a : contFrac (1/b) where
>       (a,b) = properFraction x
> 
>     -- turn truncated continued fraction into a rational number
>     unFrac :: Integral a => [a] -> Ratio Integer
>     unFrac [c] = fromIntegral c
>     unFrac (c:cs) = fromIntegral c + 1/unFrac cs
> 
>     -- make rational from first n terms of cont frac of a real
>     approx :: RealFrac a => Int -> a -> Ratio Integer
>     approx n x = unFrac . take n . contFrac $ x
> 
>     -- print first 4 rational approximations of r
>     main = forM_ [1..4] $ \i -> do
>               print (approx i r)
> 
> The output for the first 4 terms is:
> 
>     *Main> main
>     0 % 1
>     1 % 13
>     8 % 105
>     70368744177665 % 923589767331853
> 
> You can see the 4th approximant is basically noise.  The 5th is much
> worse noise, several dozen digits.

Ok. Here, you should truncate the search when the fraction gives the same number of significant digits in r, as input to the problem, which should occur at 8/105.

7.61904761904762e-02 8e 105e f/ f- fs.
0.00000000000000e+00  ok

Krishna

[toc] | [prev] | [next] | [standalone]

#57601

On Wednesday, April 19, 2017 at 12:41:10 AM UTC-7, Paul Rubin wrote:
> If you want to find the ratio by continued fractions, here's Haskell
> code that I just wrote:

Way back in 1984 when I was learning Forth, I was very impressed by */ that I read about in: "Starting Forth." Of course, I wondered how the ratios were obtained. What if you need a ratio that isn't listed in the book's short list of common ratios? Most likely, Charles Moore wrote a program in Fortran to calculate the ratios. That is a sad commentary on Forth that the inventor of Forth has to switch to Fortran to obtain the ratios that he uses in Forth with his */ that is arguably one of the coolest features of Forth (and which is lacking in Fortran). Now I see that you have written a program in Haskell to calculate the ratios. All the Forth enthusiasts jump to another language when they encounter flaws in Forth.

Later that year, after my high-school graduation, an article appeared in "Forth Dimensions" that described how to write continued fractions in Forth. This code required significant amounts of support, including writing UD/MOD in high-level Forth. At the time, everybody knew that Forth-83 had serious problems due to the fact that Charles Moore had been kicked out of Forth Inc. in 1982 and hence Forth-83 was written by Elizabeth Rather. The expectation however, was that a new standard would be coming soon --- a standard designed by Forth programmers, rather than a salesperson. A lot of time passed, and then ANS-Forth came out in 1994, but it was no better than Forth-83. Later on, in the 21st century, I upgraded the 1984 program to run under ANS-Forth --- but ANS-Forth still has the same problems as Forth-83 --- the only difference is that ANS-Forth is available on 32-bit computers.

[toc] | [prev] | [next] | [standalone]

#57655

In article <718428df-9808-4765-a9a1-48d43a971d24@googlegroups.com>,
 <krishna.myneni@ccreweb.org> wrote:
>On Tuesday, April 18, 2017 at 9:40:04 PM UTC-5, Paul Rubin wrote:
>> krishna.myneni@ccreweb.org writes:
>> > Rosetta Code has a number of entries for code to converting a
>> > (computer representation) of a real number into a rational number, or
>> > to the nearest expressible fraction, a/b, where both a and b are
>> > integers. There is even a Forth entry!
>>=20
>> The Forth entry looks pretty naive unfortunately. ...
>
>Most of the entries in other languages look naive also, in the sense that t=
>hey appear to be brute force searches.
>
>> I'm not sure I
>> totally understood the problem but I thought the sequence of
>> lowest-terms rational approximations to a real was given by the
>> continued fraction sequence.
>>=20
>> > r =3D -sqrt(2/(3^3*7)).
>> > When the finite representation of the real value is given, say to
>> > double precision, a smart algorithm to find a and b, or the prime
>> > powers representation would be handy.=20
>>=20
>> If for some reason you didn't want to do the calculation directly, like
>> you're on an MCU or something, maybe you could keep a table of the
>> logarithms of a list of primes.  Then you can find -log(r) by adding
>> suitable multiples of the needed primes, and dividing by 2.
>
>Just so the problem is clear, I have a double precision value for r, which =
>means I know it to about 16 significant digits. I know r, with infinite pre=
>cision, is a rational number, but my machine representation is a rounded ve=
>rsion to the limited precision. Now, I want to find "a" and "b" where r =3D=
> a/b. Dividing the floating point representations of a and b, ra and rb, sh=
>ould give me a number equal to r in the machine representation.
>


Let's see ...

~/PROJECT/ciforth$ lina64 -e
AMDX86 ciforth beta 2017Apr18
...
COPYRIGHT (c) ...
This program is free ... redistribute .. blabla.

 WANT -fp-
...
OK
FINIT
OK
HEX PI
OK
FDUP 100 F.R
3.243F6A8885A308D00000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000
000000000000000000000000000000000000000000000000000000000_0  OK
\ I meant 100 not 256 ...
3.243F6A8885A308D .S

S[  3243F6A8885A308D 0 ] OK
\ Fits in a single.
----------
So apparently
3243F6A8885A308D : 1000000000000000
is the ratio (for this particular fp : Intel 8087)

Where is the beef?

The refusal to use hex for floats make them seem sooooo mysterious.
(Insisting on E for exponent while in 1968 we knew better.)

\ Continuing:
3243F6A8885A308D  1000000000000000 DECIMAL
 OK
 BEGIN OVER /MOD . DUP WHILE SWAP REPEAT
 0 3 7 15 1 292 1 1 1 2 1 3 1 14 2 1 1 2 2 1 1 2 3 1 1 6 1 1 5 7
1 1 23 2 6 1 5 3 1 2  OK
.
0 OK

There you have your coefficients, the corresponding continued fraction
will represent the computer fp exactly, because it is a rational.

(Unless you're clairvoyant and *know* that number is supposed to be
the mathematical pi, that is all their is to it.
If you know it is pi, you can find an equally good approximation of pi
with half the coefficients.
You cannot find a better approximation for
3243F6A8885A308D / 1000000000000000
because it is exact. )

>Krishna
> =20
-- 
Albert van der Horst, UTRECHT,THE NETHERLANDS
Economic growth -- being exponential -- ultimately falters.
albert@spe&ar&c.xs4all.nl &=n http://home.hccnet.nl/a.w.m.van.der.horst

[toc] | [prev] | [next] | [standalone]

#57672

On Thursday, April 20, 2017 at 1:14:23 PM UTC-5, Albert van der Horst wrote:
...
> 
> S[  3243F6A8885A308D 0 ] OK
> \ Fits in a single.
> ----------
> So apparently
> 3243F6A8885A308D : 1000000000000000
> is the ratio (for this particular fp : Intel 8087)
> 

What does "S[  3243F6A8885A308D 0 ]" do?

> ...
> \ Continuing:
> 3243F6A8885A308D  1000000000000000 DECIMAL
>  OK
>  BEGIN OVER /MOD . DUP WHILE SWAP REPEAT
>  0 3 7 15 1 292 1 1 1 2 1 3 1 14 2 1 1 2 2 1 1 2 3 1 1 6 1 1 5 7
> 1 1 23 2 6 1 5 3 1 2  OK
> .
> 0 OK
> 
> There you have your coefficients, the corresponding continued fraction
> will represent the computer fp exactly, because it is a rational.
> 
> (Unless you're clairvoyant and *know* that number is supposed to be
> the mathematical pi, that is all their is to it.
> If you know it is pi, you can find an equally good approximation of pi
> with half the coefficients.
> You cannot find a better approximation for
> 3243F6A8885A308D / 1000000000000000
> because it is exact. )

It is "exact" in the sense of the finite precision machine representation. However, are you sure there is no other ratio, using smaller integers, which will you give you the "exact" real number to the same precision?

In the previous example I gave,

r = 7.61904761904762e-02 

one can certainly write it as the ratio of integers, 761904781904762 : 100000000000000000 . However, the smallest pair of integers which will yield the same precision (as written above) is 8 : 105. The continued fraction can be truncated at the point where the rational approximation no longer provides any further improvement in error.

In some instances, such as for the squares of the Wigner 3-j symbols, we know in advance that the results are rational numbers. However, on a computer, they are often calculated in such a way to give a floating point number and it is difficult to find the corresponding rational number by hand. 

Krishna

[toc] | [prev] | [next] | [standalone]

#57679

krishna.myneni@ccreweb.org writes:
> In some instances, such as for the squares of the Wigner 3-j symbols,
> we know in advance that the results are rational numbers. However, on
> a computer, they are often calculated in such a way to give a floating
> point number and it is difficult to find the corresponding rational
> number by hand.

You could use mpfr and do the calculation in exact rational arithmetic
from the beginning.

[toc] | [prev] | [next] | [standalone]

#57575

El martes, 18 de abril de 2017, 8:49:31 (UTC-3), krishna...@ccreweb.org  escribió:
> Rosetta Code has a number of entries for code to converting a (computer representation) of a real number into a rational number, or to the nearest expressible fraction, a/b, where both a and b are integers. There is even a Forth entry!
> 
> https://rosettacode.org/wiki/Convert_decimal_number_to_rational
> 
> However, it seems like many of the algorithms expressed there are brute force searches, and don't make use of the fundamental theorem of arithmetic to find the suitable integers a and b (though the code which uses library functions may do so).
> 
> This particular problem has been in the back of my mind for some time, and it is one that comes up from time to time. For example, in the quantum physics of angular momentum, the Wigner 3-j symbols are real numbers between zero and one, expressible as +/-sqrt(a/b), and if one computes the quantities as real numbers, it would be nice to find (a/b) in an efficient manner. In tables of these symbols, the rational number a/b is expressed as a series of powers of prime numbers and a sign term, e.g. -1 [1 -3 0 -1] represents 
> 
> r = -sqrt(2/(3^3*7)).
> 
> When the finite representation of the real value is given, say to double precision, a smart algorithm to find a and b, or the prime powers representation would be handy. The real to rational code is certainly useful, but I think the algorithms should be able to take advantage of the fact that |r| <= 1. It's an interesting inversion problem.
> 
> Incidentally, here is Forth code to compute Wigner 3j symbols -- it is somewhat limited in the range of usable angular momentum quantum numbers, though.
> 
> ftp://ccreweb.org/software/kforth/examples/fsl/extras/cg.4th
> 
> and accompanying test code,
> 
> ftp://ccreweb.org/software/kforth/examples/fsl/extras/cg-test.4th
> 
> Krishna

I convert from fixed point 16.16 to fraction for a calculator

the complete lib in
https://github.com/phreda4/reda4/blob/master/r4/Lib/frac.txt

the lib has 50 lines, make a fraction of decimal part symplify and sum de integer par, I guess.



::gcd | n1 n2 -- gcd
	( 1? )(	swap over mod ) drop ;

::fsimp | n1 n2 -- n1 n2
	2dup gcd 0? ( drop ; ) 
	>r r / swap r> / swap
	-? ( swap neg swap neg ) ;

::fix2f | fix -- n1 d1
	dup 16 >> swap $10000 fsimp rot over * rot + swap ;

[toc] | [prev] | [next] | [standalone]

#57577

On Tuesday, April 18, 2017 at 9:58:24 PM UTC-5, Pablo Hugo Reda wrote:
...
> I convert from fixed point 16.16 to fraction for a calculator
> 
> the complete lib in
> https://github.com/phreda4/reda4/blob/master/r4/Lib/frac.txt
> 
> the lib has 50 lines, make a fraction of decimal part symplify and sum de integer par, I guess.
> 
> 
> 
> ::gcd | n1 n2 -- gcd
> 	( 1? )(	swap over mod ) drop ;
> 
> ::fsimp | n1 n2 -- n1 n2
> 	2dup gcd 0? ( drop ; ) 
> 	>r r / swap r> / swap
> 	-? ( swap neg swap neg ) ;
> 
> ::fix2f | fix -- n1 d1
> 	dup 16 >> swap $10000 fsimp rot over * rot + swap ;

Hi Pablo,

Can you give me an example of using this to solve the originally stated problem?

Krishna

[toc] | [prev] | [next] | [standalone]

#57596

El miércoles, 19 de abril de 2017, 0:05:41 (UTC-3), krishna...@ccreweb.org  escribió:
> On Tuesday, April 18, 2017 at 9:58:24 PM UTC-5, Pablo Hugo Reda wrote:
> ...
> > I convert from fixed point 16.16 to fraction for a calculator
> > 
> > the complete lib in
> > https://github.com/phreda4/reda4/blob/master/r4/Lib/frac.txt
> > 
> > the lib has 50 lines, make a fraction of decimal part symplify and sum de integer par, I guess.
> > 
> > 
> > 
> > ::gcd | n1 n2 -- gcd
> > 	( 1? )(	swap over mod ) drop ;
> > 
> > ::fsimp | n1 n2 -- n1 n2
> > 	2dup gcd 0? ( drop ; ) 
> > 	>r r / swap r> / swap
> > 	-? ( swap neg swap neg ) ;
> > 
> > ::fix2f | fix -- n1 d1
> > 	dup 16 >> swap $10000 fsimp rot over * rot + swap ;
> 
> Hi Pablo,
> 
> Can you give me an example of using this to solve the originally stated problem?
> 
> Krishna

Krishna

I found a bug!,not work for number with not zero integer part, the code for test is

:main
	show clrscr

		0.6190 | 7.619 not work
		dup "%f " print  | print fixpoint
		fix2f 2dup       | convert to fraction
		swap "%d/%d" print | print in fraction
    	        16 <</ " %f" print | back to fixed point

		'exit >esc<
		cminiflecha ;

I don't remenber how code this.
for 0.619 work ok but for 7.619 give me 14.619, the integer part is doubled! 
sorry.

[toc] | [prev] | [next] | [standalone]

Page 2 of 4 — ← Prev page 1 [2] 3 4  Next page →

Back to top | Article view | comp.lang.forth

csiph-web