|
===Plan===
This is F77 style, and directly translates the given calculation as per ''formula translation''. It turns out that the normal integers suffice for the demonstration, except for just one of the products of the three primes: 41x1721x35281 = 2489462641, which is bigger than 2147483647, the 32-bit limit. Fortunately, INTEGER*8 variables are also available, so the extension is easy. Otherwise, one would have to mess about with using two integers in a bignum style, one holding say the millions, and the second the number up to a million.
===Integer arithmetic===
Integer arithmetic can be a delicate matter, not just because of overflows. With division, sometimes the formula will never produce a remainder, as in n(n + 1)/2, and sometimes the divisor is always a factor for more complex reasons - as with <code>(P1 - 1)*(H3 + P1)/D</code>. But if a term produces a fractional part, is it to be discarded or not? Some languages (such as Pascal) insist on "div" instead of "/" so that one thinks about truncating division, others (such as pl/i) just generate the fractional part and proceed: perhaps the result will round? truncate? to the desired value, or perhaps other terms will cancel out the fractional parts and all will be well, or, perhaps they will combine and the result will be out by one, ''etc''. And some formulae require the fractional parts, as in the expression for the n'th Fibonacci number where each term generates an ''infinite'' non-repeating fractional part, and, they cancel ''exactly'', even for large values of n... In source code,
:<code>(((1 + SQRT(5))/2)**N - ((1 - SQRT(5))/2)**N)/SQRT(5)</code>
or, in mathematical notation (possibly not rendered),
:<math>\frac{1}{\sqrt{5}}\left[\left(\frac{1 + \sqrt{5}}{2}\right)^n - \left(\frac{1 - \sqrt{5}}{2}\right)^n\right]</math>
===What model MOD is your MOD?===
Directly translating the given formulae involves a trap. There is an expression ''-Prime1 squared mod h3'', and translating this to <code>MOD(-P1**2,H3)</code> may not work. The behaviour of the MOD function for negative numbers varies between language, compiler, and cpu. For positive numbers there is no confusion. However, in some cases one wants an affine shift as in day-of-the-week calculations from a daynumber, D: given MOD(D,7), one wants the same result with MOD(D - 7,7) or any other such shift and indeed, this happens with some versions of MOD. But not all, as Prof. Knuth remarks in his description of the calculation of the date for Easter. The MOD function can also work as follows:
<pre>
D -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7
MOD(D,3) -1 0 -2 -1 0 -2 -1 0 1 2 0 1 2 0 1
MOD(3 + MOD(D,3),3) 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1
</pre>
Thus, <code>MOD(-1,3) + MOD(4,3) = -1 + 1</code> and <code>MOD(3,3) = 0</code> so that ''(a + b) mod n'' does come out as ''[(a mod n) + (b mod n)] mod n'' as is expected, but in a different way: (-1 ''mod'' 3) + (4 ''mod'' 3) = 2 + 1, thus the need for the second ''mod''. So, MOD and ''mod'' are different. However, this MOD is consistent with integer division that truncates. If ''q·f = n/d'' where q is the integer part and f the fractional part of the division, then the remainder is ''fd'' and has a sign too. Another computer/compiler/language version of MOD may not differ from ''mod''.
Here I am supposing that ''(h3+Prime1)*(Prime1-1) mod d == 0 and -Prime1 squared mod h3 == d mod h3'' is evaluated as ''{[(h3 + Prime1)*(Prime1 - 1) mod d] == 0} and {[-(Prime1 squared) mod h3] == d mod h3}'' which is to say it is ''[-(Prime1 squared)] mod h3'' not ''-[(Prime1 squared) mod h3]'' where the first case involves the MOD(negative number) while the second is -MOD(positive number) and the only way that result could match ''d mod h3'' is when both are zero.
Testing with the peculiar sign ignored via <code>.AND. (MOD(P1**2,H3) .EQ. MOD(D,H3))) THEN</code> gives limited results:
<pre>
Carmichael numbers that are the product of three primes:
P1 x P2 x P3 = C
3 11 17 561
5 29 73 10585
7 19 67 8911
13 61 397 314821
13 37 241 115921
19 43 409 334153
31 991 15361 471905281
37 109 2017 8134561
41 1721 35281 2489462641
43 631 13567 368113411
43 271 5827 67902031
43 127 2731 14913991
61 421 12841 329769721
61 181 5521 60957361</pre>
Which turns out to be when the two terms have remainders of one.
===Source===
So, using the double MOD approach - which gives the same result for either style of MOD... <lang Fortran> LOGICAL FUNCTION ISPRIME(N) !Ad-hoc, since N is not going to be big...
| 