-4

This is math calculation that i need to program - define sequence an as

  • a1 = 1
  • a(n+1) = 1 / (2 × [an] - an + 1) where [x] is the floor function

If a2014 = p/q, find p + q

And, this is the code I tried:

#include <stdio.h>
#include <math.h>

main() {
   int n=1;
   double an=1.0;
   double an_floor=0.0;

   while (n<=2014) {
      an_floor = floor(an);
      an = 1 / (2*an_floor-an+1);
      n = n + 1;
   }

   printf("%lf", an);

   return 0;
}

Problems;

  • I cannot compile (using web complier such as codepad, etc.)

  • Don't know how to make result

  • Fraction result is wrong

Sinan Ünür
  • 116,958
  • 15
  • 196
  • 339
John Jang
  • 17

2 Answers2

1

For someone who always likes to remind people about floating point errors, I was caught asleep at the wheel in this case.

If you run the following program,

/* a1 = 1, a(n+1) = 1/(2*[an]-an+1) ([x] is floor function)
 * a2014 = p/q
 * find p+q
 */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#define LAST_ELEMENT 2014

static double
next_element(double prev) {
    return 1.0 / (2.0 * floor(prev) - prev + 1.0);
}


int main(void) {
    int i = 0;
    double last = 1.0;

    for (i = 0; i < LAST_ELEMENT; i += 1) {
        last = next_element( last );
    }

    printf("%g\n", last);
    return EXIT_SUCCESS;
}

you get this output:

C:\...\Temp> cl /fp:precise /O2 rrr.c
C:\...\Temp> rrr.exe
2.25

But, that is due to floating point error as @JJoao points out. He outlines a specific way you can deal with this particular issue.

Another way is to utilize an arbitrary precision math library to help you. First, let's verify the issue using a quick Perl script:

#!/usr/bin/env perl

use strict;
use warnings;

use POSIX qw( floor );

sub next_element { '1' / (('2' * floor($_[0])) - $_[0] + '1.0') }

sub main {
    my ($x, $n) = @_;
    $x = next_element( $x ) for 1 .. $n;
    printf("%g\n", $x);
}

main(1, 2014);

Output:

C:\...\Temp> perl t.pl
2.25

Same as the C program.

Now, use Perl's bignum:

C:\...\Temp> perl -Mbignum t.pl
5.83333

This increased accuracy comes at a performance cost. Without bignum, the script runs in 0.125 seconds with about 0.094 is spent doing other things than the calculation. With bignum, it takes about two seconds.

Now, modern C provides various facilities for manipulating how floating numbers are rounded. In this particular case, given the nature of the floor function, setting the rounding mode to FE_DOWNWARD will fix the issue:

#include <fenv.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#define LAST_ELEMENT 2014

static double
next_element(double prev) {
    return 1.0 / (2.0 * floor(prev) - prev + 1.0);
}


int main(void) {
    int i = 0;
    double last = 1.0;

    const int original_rounding = fegetround();
    fesetround(FE_DOWNWARD);

    for (i = 0; i < LAST_ELEMENT; i += 1) {
        last = next_element(last);
    }

    fesetround(original_rounding);
    printf("%g\n", last);

    return EXIT_SUCCESS;
}

Output:

C:\...\Temp> cl /O2 /fp:precise rrr.c            
/out:rrr.exe                                                           

C:\...\Temp> rrr                                 
5.83333
Community
  • 1
  • 1
Sinan Ünür
  • 116,958
  • 15
  • 196
  • 339
  • would it be `2.2518e+15` = infinty? – JJoao May 07 '15 at 11:45
  • **Very interesting!!!** 1) in unix I can not reproduce the behavior of the C version. (I get `2.5` (I dont know what the correspondent to `/fp:precise` option) – JJoao May 08 '15 at 16:52
  • Hmmmm ... can you try adding `#pragma STDC FENV_ACCESS ON`? – Sinan Ünür May 08 '15 at 16:55
  • Interesting. I'll try to figure it for `gcc` when I get a chance. – Sinan Ünür May 08 '15 at 17:01
  • 1
    replace `FE_DOWNWARD` by `FE_UPWARD` and ti works! It is curious that we already lost 5 decimal digits in the result precision... 5.8333333333428480 5.83333 – JJoao May 08 '15 at 17:21
  • 1
    Sinan, it looks like "excess floating-point precision", "rounding" and similar are a difficult areas, (bugs, possible differences bw compilers, etc). In our problem: (1) `fesetround(FE_UPWARD)` fixes the floor problem (Excellent) - and probably with a minor increase of the global final error (very minor problem). (2) `perl -Mbignum` works, and is a very clear and generic technique (Thank you for telling me about big* perl modules). (3) `perl -Mbigrat` works and is a no representation error solution. – JJoao May 09 '15 at 12:13
1

This is a dangerous sequence:

floor(value with potential representation error) = Big errors!!

floor(3-epsilon)= 2
floor(3+epsilon)= 3

Big number of cumulative real number calculation = Big errors!!

You need to study the numeric-precision errors!

using the OP code we get

1
0.5
2
0.333333
1.5
0.666667
3
0.5           <===== Big precision ERROR it should be 0.25
2
1
2.2518e+15A   <===== infinity
4.44089e-16   <===== zero
1
0.5
....(wrong) cyclic behavior

final (wrong!!) result = 2.2518e+15A <===== infinity

Update: But what should we do?

1) floor problem: substitute floor(x) by floor(x + 2*expectedError)

2) about cumulative errors: Replace double representation by a (v1:int, v2:int) fraction representation.An -> v1/v2

newAn = newv1/newv2 =
      = 1 / (2 * flo(An) + 1 - An) =
      = 1 / (2 * flo(v1/v2) + 1 - v1/v2) =
      = ... =
      = v2 / ( 2 * v2 * flo(v1/v2) + v2 - v1)

This way

newv1 = v2
newv2 =  2 * v2 * flo(v1/v2) + v2 - v1

In C:

main() {
   int v1=1, v2=1, n=1, oldv1;
   double an_floor;

   while (n<=2014) {
     an_floor = floor((0.0+v1)/v2 + 0.0000000001);

     oldv1=v1;
     v1=v2;
     v2=an_floor * v2 * 2 + v2 - oldv1;
     n++;
     //  printf("%g\n", (0.0+v1)/v2);
  }

  printf("%d/%d=%g\n", v1, v2, (0.0+v1)/v2);   //      35/6=5.83333
}

\thanks{Sinan Ünür}

JJoao
  • 4,891
  • 1
  • 18
  • 20