1

The problem I'm attempting to tackle at the moment involves computing the order of 10 modulo(n), where n could be any number less than 1000. I have a function to do exactly that, however, I am unable to obtain accurate results as the value of the order increases.

The function works correctly as long as the order is sufficiently small, but returns incorrect values for large orders. So I stuck in some output to the terminal to locate the problem, and discovered that when I use exponentiation, the accuracy of my reals is being compromised.

I declared ALL variables in the function and in the program I tested it from as real(kind=nkind) where nkind = selected_real_kind(p=18, r=308). Any numbers explicitly referenced are also declared as, for example, 1.0_nkind. However, when I print out 10**n for n counting up from 1, I find that at 10**27, the value is correct. However, 10**28 gives 9999999999999999999731564544. All higher powers are similarly distorted, and this inaccuracy is the source of my problem.

So, my question is, is there a way to work around the error? I don't know of any way to use a more extended precision than I'm already using in the calculations.

Thanks, Sean

*EDIT: There's not much to see in the code, but here you go:

integer, parameter :: nkind = selected_real_kind(p=18, r = 308)
real(kind=nkind) function order_ten_modulo(n)

real(kind=nkind) :: n, power
power = 1.0_nkind

if (mod(n, 5.0_nkind) == 0 .or. mod(n, 2.0_nkind) == 0) then
    order_ten_modulo = 0
    return
end if


do 
    if (power>300.0) then ! Just picked this number as a safeguard against endless looping -
        exit
    end if

    if (mod(10.0_nkind**power, n) == 1.0_nkind) then
    order_ten_modulo = power
    exit
    end if
    power = power + 1.0_nkind
end do

return

end function order_ten_modulo
Sean McBane
  • 21
  • 3
  • Project Euler #26 deals with recurring digits of `1/n`, which I have a hard time connecting to your issue. There is no way of reducing imprecision of `1**27` in Fortran; but there might be a way to not require such a calculation, doing things another way. Finally, even for well-known sites like Project Euler, please note the Stack Overflow policy - your question should be self-contained, understandable without needing to open another page. (Linking for convenience, such as including code *and* linking to online demo, is welcome and encouraged.) – Amadan Apr 17 '14 at 04:56
  • Thanks for the thought - the issue stems from what I found while investigating the issue: the period of the repeating fraction for a prime number n (in base 10) is equal to the order of 10 modulo n. And by playing around with it a little bit, I found that rules involving the prime factors can be used for numbers that aren't prime... Anyway, since that order is the smallest k such that 10^k = 1 (modulo n), I just wrote a function that looked for that k. But I'll investigate a more elegant algorithm. – Sean McBane Apr 17 '14 at 05:12
  • This might be most elegant mathematically. But algorithmically, with data structures that you have available, I would suggest you try mimicking the elementary school long division, step by step. This way, you only need integers from `0` to `1000`, and no floats at all. – Amadan Apr 17 '14 at 05:23
  • 2
    Most Fortran compilers offer quadrupole precision. Or use an arbitrary precision library: http://stackoverflow.com/questions/21792114/arbitrary-precision-arthimetic-in-fortran – M. S. B. Apr 17 '14 at 07:52
  • instead of describing the program, please post some actual code. – steabert Apr 17 '14 at 09:10
  • @steabert: I added the function in question to the original post. It's part of a module with various other mathematical routines; the parameter declaration is at the beginning of the module, of course. – Sean McBane Apr 17 '14 at 14:24
  • @Amadan: I thought of that, since I've done something similar for working with numbers hundreds of digits long before, but I can't think of a very simple way to recognize the length of the repeating pattern, when it may be hundreds of digits long. – Sean McBane Apr 17 '14 at 14:34
  • @M.S.B.: Thanks! Those libraries are just the sort of thing I was looking for. The quadruple precision type is accurate up to 10^48, which would get me considerably further, but I think I still probably need the capabilities from a library. Thanks again. – Sean McBane Apr 17 '14 at 14:38
  • 1
    @SeanMcBane: I am not a maths expert, but I believe the pattern will start to repeat when, in long division, you get the same number (modulo n) as a dividend. Thus, an array of length `n` suffices to memorise what you've seen before. – Amadan Apr 17 '14 at 23:43

0 Answers0