the secret life of the mod function

define a function for real numbers x & y, y != 0 as follows:

mod(x, y) = x - floor(x/y)*y = frac(x/y)*y

so, when x = 10 and y = 3, we have mod(10,3) = 10 - 3*3 = 1 and when x = 10 and y = pi, mod(10,pi) = 10-3*pi.

everything is fine mathematically, but suppose we want to compute the mod function for x & y floating point numbers accurate to 30 significant figures. suppose, also the floating point numbers can have exponents up to 1e100.

define another function,    modtwopi(x) = mod(x, 2*pi)

there are immediately two problems with the modtwopi function, where a 30 digit x value can be supplied and an accurate 30 digit answer is expected.

bulletfirstly, somewhere internally, 2pi needs to be stored to 100 decimal places.

consider modtwopi(1e50) = frac(1e50/2pi)*2pi =
frac(15915494309189533576888376337251436203445964574045.644874766734405889679763422653... )*2pi
= 0.644874766734405889679763422653*2pi
= 4.051867659316482244826000541880
(or -2.231317647863104232099286224678 as given by scicalc because it adjusts to -pi to +pi)

notice how, in this example, the 50 digit integer part is discarded and 30 more correct digits are required after the decimal point. these are multiplied up to give the final answer. here 80 digits of 1/2pi are needed.

bulletsecondly, as you can see from the example above, double the precision arithmetic is required for intermediate calculation. if we want 30 digits we will need 60 temporary digits.

the mod function itself suffers from the double precision problem. the modtwopi function has the additional trouble of having 2pi implied as the divisor and therefore requires a high precision 2pi representation.

ive left some gaps because the above example would suggest that more than 100 digits of 1/2pi are needed and more than 60 digits of arithmetic are needed. let look at this now.

because we are working to 30 (and in decimal), we can represent x as n*10^b where n is an integer and b = log10(x)-30  (the exponent minus 30). and mod(x, y) = mod(n*mod(10^b),y) breaking the problem into two parts. we need 30 digit precision for mod(10^b,y) but using an auxiliary table of the digits of 1/2pi, then 60 for the n*mod(10^b,y) and 60 for the final mod, yielding a 30 digit remainder from the division.

note that mod(a*b, y) != mod(mod(a, y) * mod(b, y), y) unless both a & b are integers.

what use is this?

many calculators no longer reject large inputs to trig functions. take sin(x), for example. what is sin(1e50) (radians)?

sin(1e50) = sin(modtwopi(1e50)) = sin(-2.231317647863104232099286224678) = -.7896724934293100827102895399174 correct to 30 places. either they should reject the input or get it right. answers from various models are tabulated here. even for moderate values of x = 10,000. calculators will have trouble getting all their digits correct since 4 extra places of pi are required internally (in this case).

lets test this using two machines, the commodore SR414R, a straightforward scientific and the HP 32sii

calculator sin(1.0) radians sin(10,000) radians
reference value 0.841470984807897 -0.30561438888825
commodore SR4148R 0.841470984, all 9 digits correct -0.305614358, only 7 correct now.
HP 32Sii 0.841470984808, all 12 digits correct -0.30561438888, all 12 digits correct

the commodore gets two digits less than for values within the normal trig domain of -pi to +pi. the hp, however, does remarkably well. perhaps it has an extra precision mod twopi.