Usages:

Right click on identifier in code to find usages!

 
1   package org.apfloat;
2   
3   /**
4    * 
Various mathematical functions for arbitrary precision integers.
5    *
6    * 
@version 1.6
7    * 
@author Mikko Tommila
8    */

9   
10  public class ApintMath
11  {
12      private ApintMath()
13      {
14      }
15  
16      /**
17       * 
Integer power.
18       *
19       * 
@param x Base of the power operator.
20       * 
@param n Exponent of the power operator.
21       *
22       * 
@return <code>x</code> to the <code>n</code>:th powerthat is <code>x<sup>n</sup></code>.
23       *
24       * 
@exception java.lang.ArithmeticException If both <code>x</code> and <code>n</code> are zero.
25       */

26  
27      public static Apint pow(Apint xlong n)
28          throws ArithmeticExceptionApfloatRuntimeException
29      {
30          if (n == 0)
31          {
32              if (x.signum() == 0)
33              {
34                  throw new ArithmeticException("Zero to power zero");
35              }
36  
37              return new Apint(1, x.radix());
38          }
39          else if (n < 0)
40          {
41              return Apint.ZERO;
42          }
43  
44          // Algorithm improvements by Bernd Kellner
45  
        int b2pow = 0;
46  
47          while ((n & 1) == 0)
48          {
49              b2pow++;
50              n >>= 1;
51          }
52  
53          Apint r = x;
54  
55          while ((n >>= 1) > 0)
56          {
57              x = x.multiply(x);
58              if ((n & 1) != 0)
59              {
60                  r = r.multiply(x);
61              }
62          }
63  
64          while (b2pow-- > 0)
65          {
66              r = r.multiply(r);
67          }
68  
69          return r;
70      }
71  
72      /**
73       * 
Square root and remainder.
74       *
75       * 
@param x The argument.
76       *
77       * 
@return An array of two apints<code>[qr]</code>where <code>q<sup>2</sup> + r = x</code>.
78       *
79       * 
@exception java.lang.ArithmeticException If <code>x</code> is negative.
80       */

81  
82      public static Apint[] sqrt(Apint x)
83          throws ArithmeticExceptionApfloatRuntimeException
84      {
85          return root(x, 2);
86      }
87  
88      /**
89       * 
Cube root and remainder.
90       *
91       * 
@param x The argument.
92       *
93       * 
@return An array of two apints<code>[qr]</code>where <code>q<sup>3</sup> + r = x</code>.
94       */

95  
96      public static Apint[] cbrt(Apint x)
97          throws ApfloatRuntimeException
98      {
99          return root(x, 3);
100     }
101 
102     /**
103      * 
Positive integer root and remainder.<p>
104      *
105      * 
Returns the <code>n</code>:th root of <code>x</code>,
106      * 
that is <code>x<sup>1/n</sup></code>rounded towards zero.
107      *
108      * 
@param x The argument.
109      * 
@param n Which root to take.
110      *
111      * 
@return An array of two apints<code>[qr]</code>where <code>q<sup>n</sup> + r = x</code>.
112      *
113      * 
@exception java.lang.ArithmeticException If <code>n</code> and <code>x</code> are zeroor <code>x</code> is negative and <code>n</code> is even.
114      */

115 
116     public static Apint[] root(Apint xlong n)
117         throws ArithmeticExceptionApfloatRuntimeException
118     {
119         if (n == 0)
120         {
121             if (x.signum() == 0)
122             {
123                 throw new ArithmeticException("Zeroth root of zero");
124             }
125 
126             Apint one = new Apint(1, x.radix());
127             return new Apint[] { onex.subtract(one) };
128         }
129         else if (x.signum() == 0)
130         {
131             return new Apint[] { xx };                        // Avoid division by zero
132 
        }
133         else if (x.equals(Apint.ONE) || n == 1)
134         {
135             return new Apint[] { xApint.ZERO };
136         }
137         else if (n < 0)
138         {
139             return new Apint[] { Apint.ZEROx };               // 1 / x where x > 1
140 
        }
141 
142         long precision = x.scale() / n + Apint.EXTRA_PRECISION;
143         Apfloat approxX = x.precision(precision);
144         Apfloat approxRoot;
145 
146         approxRoot = ApfloatMath.root(approxXn);
147 
148         Apint root = approxRoot.truncate(),                             // May be one too big or one too small
149 
              pow = pow(rootn);
150 
151         if (abs(pow).compareTo(abs(x)) > 0)
152         {
153             // Approximate root was one too big
154 

155             pow = (x.signum() >= 0 ? powXMinus1(powrootn) : powXPlus1(powrootn));
156             root = root.subtract(new Apint(x.signum(), x.radix()));
157         }
158         else
159         {
160             // Approximate root was correct or one too small
161 

162             Apint powPlus1 = (x.signum() >= 0 ? powXPlus1(powrootn) : powXMinus1(powrootn));
163 
164             if (abs(powPlus1).compareTo(abs(x)) <= 0)
165             {
166                 // Approximate root was one too small
167 

168                 pow = powPlus1;
169                 root = root.add(new Apint(x.signum(), x.radix()));
170             }
171         }
172 
173         Apint remainder = x.subtract(pow);
174 
175         assert (remainder.signum() * x.signum() >= 0);
176 
177         return new Apint[] { rootremainder };
178     }
179 
180     private static Apint powXMinus1(Apint powApint xlong n)
181         throws ApfloatRuntimeException
182     {
183         Apint one = new Apint(1, x.radix());
184 
185         if (n == 2)
186         {
187             // (x - 1)^2 = x^2 - 2*x + 1
188 
            pow = pow.subtract(x).subtract(x).add(one);
189         }
190         else if (n == 3)
191         {
192             // (x - 1)^3 = x^3 - 3*x^2 + 3*x - 1 = x^3 - 3*x*(x - 1) - 1
193 
            pow = pow.subtract(new Apint(3, x.radix()).multiply(x).multiply(x.subtract(one))).subtract(one);
194         }
195         else
196         {
197             pow = pow(x.subtract(one), n);
198         }
199 
200         return pow;
201     }
202 
203     private static Apint powXPlus1(Apint powApint xlong n)
204         throws ApfloatRuntimeException
205     {
206         Apint one = new Apint(1, x.radix());
207 
208         if (n == 2)
209         {
210             // (x + 1)^2 = x^2 + 2*x + 1
211 
            pow = pow.add(x).add(x).add(one);
212         }
213         else if (n == 3)
214         {
215             // (x + 1)^3 = x^3 + 3*x^2 + 3*x + 1 = x^3 + 3*x*(x + 1) + 1
216 
            pow = pow.add(new Apint(3, x.radix()).multiply(x).multiply(x.add(one))).add(one);
217         }
218         else
219         {
220             pow = pow(x.add(one), n);
221         }
222 
223         return pow;
224     }
225 
226     /**
227      * 
Returns an apint whose value is <code>-x</code>.
228      *
229      * 
@deprecated Use {@link Apint#negate()}.
230      *
231      * 
@param x The argument.
232      *
233      * 
@return <code>-x</code>.
234      */

235 
236     @Deprecated
237     public static Apint negate(Apint x)
238         throws ApfloatRuntimeException
239     {
240         return x.negate();
241     }
242 
243     /**
244      * 
Absolute value.
245      *
246      * 
@param x The argument.
247      *
248      * 
@return Absolute value of <code>x</code>.
249      */

250 
251     public static Apint abs(Apint x)
252         throws ApfloatRuntimeException
253     {
254         if (x.signum() >= 0)
255         {
256             return x;
257         }
258         else
259         {
260             return x.negate();
261         }
262     }
263 
264     /**
265      * 
Copy sign from one argument to another.
266      *
267      * 
@param x The value whose sign is to be adjusted.
268      * 
@param y The value whose sign is to be used.
269      *
270      * 
@return <code>x</code> with its sign changed to match the sign of <code>y</code>.
271      *
272      * 
@since 1.1
273      */

274 
275     public static Apint copySign(Apint xApint y)
276         throws ApfloatRuntimeException
277     {
278         if (y.signum() == 0)
279         {
280             return y;
281         }
282         else if (x.signum() != y.signum())
283         {
284             return x.negate();
285         }
286         else
287         {
288             return x;
289         }
290     }
291 
292     /**
293      * 
Multiply by a power of the radix.
294      * 
Any rounding will occur towards zero.
295      *
296      * 
@param x The argument.
297      * 
@param scale The scaling factor.
298      *
299      * 
@return <code>x * x.radix()<sup>scale</sup></code>.
300      */

301 
302     public static Apint scale(Apint xlong scale)
303         throws ApfloatRuntimeException
304     {
305         return ApfloatMath.scale(xscale).truncate();
306     }
307 
308     /**
309      * 
Quotient and remainder.
310      *
311      * 
@param x The dividend.
312      * 
@param y The divisor.
313      *
314      * 
@return An array of two apints<code>[quotientremainder]</code>that is <code>[x / yx % y]</code>.
315      *
316      * 
@exception java.lang.ArithmeticException In case the divisor is zero.
317      */

318 
319     public static Apint[] div(Apint xApint y)
320         throws ArithmeticExceptionApfloatRuntimeException
321     {
322         if (y.signum() == 0)
323         {
324             throw new ArithmeticException("Division by zero");
325         }
326         else if (x.signum() == 0)
327         {
328             // 0 / x = 0
329 
            return new Apint[] { xx };
330         }
331         else if (y.equals(Apint.ONE))
332         {
333             // x / 1 = x
334 
            return new Apint[] { xApint.ZERO };
335         }
336 
337         long precision;
338         Apfloat txty;
339         Apint abqr;
340 
341         a = abs(x);
342         b = abs(y);
343 
344         if (a.compareTo(b) < 0)
345         {
346             return new Apint[] { Apint.ZEROx };  // abs(x) < abs(y)
347 
        }
348         else
349         {
350             precision = x.scale() - y.scale() + Apint.EXTRA_PRECISION;        // Some extra precision to avoid round-off errors
351 
        }
352 
353         tx = x.precision(precision);
354         ty = y.precision(precision);
355 
356         q = tx.divide(ty).truncate();           // Approximate division
357 

358         a = a.subtract(abs(q.multiply(y)));
359 
360         if (a.compareTo(b) >= 0)                // Fix division round-off error
361 
        {
362             q = q.add(new Apint(x.signum() * y.signum(), x.radix()));
363             a = a.subtract(b);
364         }
365         else if (a.signum() < 0)                // Fix division round-off error
366 
        {
367             q = q.subtract(new Apint(x.signum() * y.signum(), x.radix()));
368             a = a.add(b);
369         }
370 
371         r = copySign(ax);
372 
373         return new Apint[] { qr };
374     }
375 
376     /**
377      * 
Greatest common divisor.
378      * 
This method returns a positive number even if one of <code>a</code>
379      * 
and <code>b</code> is negative.
380      *
381      * 
@param a First argument.
382      * 
@param b Second argument.
383      *
384      * 
@return Greatest common divisor of <code>a</code> and <code>b</code>.
385      */

386 
387     public static Apint gcd(Apint aApint b)
388         throws ApfloatRuntimeException
389     {
390         return GCDHelper.gcd(ab);
391     }
392 
393     /**
394      * 
Least common multiple.
395      * 
This method returns a positive number even if one of <code>a</code>
396      * 
and <code>b</code> is negative.
397      *
398      * 
@param a First argument.
399      * 
@param b Second argument.
400      *
401      * 
@return Least common multiple of <code>a</code> and <code>b</code>.
402      */

403 
404     public static Apint lcm(Apint aApint b)
405         throws ApfloatRuntimeException
406     {
407         if (a.signum() == 0 && b.signum() == 0)
408         {
409             return Apint.ZERO;
410         }
411         else
412         {
413             return abs(a.multiply(b)).divide(gcd(ab));
414         }
415     }
416 
417     /**
418      * 
Modular multiplication. Returns <code>a * b % m</code>
419      *
420      * 
@param a First argument.
421      * 
@param b Second argument.
422      * 
@param m Modulus.
423      *
424      * 
@return <code>a * b mod m</code>
425      */

426 
427     public static Apint modMultiply(Apint aApint bApint m)
428         throws ApfloatRuntimeException
429     {
430         return a.multiply(b).mod(m);
431     }
432 
433     private static Apint modMultiply(Apint x1Apint x2Apint yApfloat inverseY)
434         throws ApfloatRuntimeException
435     {
436         Apint x = x1.multiply(x2);
437 
438         if (x.signum() == 0)
439         {
440             // 0 % x = 0
441 
            return x;
442         }
443 
444         long precision = x.scale() - y.scale() + Apfloat.EXTRA_PRECISION;       // Some extra precision to avoid round-off errors
445 
        Apint abt;
446 
447         a = abs(x);
448         b = abs(y);
449 
450         if (a.compareTo(b) < 0)
451         {
452             return x;                           // abs(x) < abs(y)
453 
        }
454 
455         t = x.multiply(inverseY.precision(precision)).truncate();               // Approximate division
456 

457         a = a.subtract(abs(t.multiply(y)));
458 
459         if (a.compareTo(b) >= 0)                // Fix division round-off error
460 
        {
461             a = a.subtract(b);
462         }
463         else if (a.signum() < 0)                // Fix division round-off error
464 
        {
465             a = a.add(b);
466         }
467 
468         t = copySign(ax);
469 
470         return t;
471     }
472 
473     /**
474      * 
Modular power.
475      *
476      * 
@param a Base.
477      * 
@param b Exponent.
478      * 
@param m Modulus.
479      *
480      * 
@return <code>a<sup>b</sup> mod m</code>
481      *
482      * 
@exception java.lang.ArithmeticException If the exponent is negative but the GCD of <code>a</code> and <code>m</code> is not 1 and the modular inverse does not exist.
483      */

484 
485     public static Apint modPow(Apint aApint bApint m)
486         throws ArithmeticExceptionApfloatRuntimeException
487     {
488         if (b.signum() == 0)
489         {
490             if (a.signum() == 0)
491             {
492                 throw new ArithmeticException("Zero to power zero");
493             }
494 
495             return new Apint(1, a.radix());
496         }
497         else if (m.signum() == 0)
498         {
499             return m;                           // By definition
500 
        }
501 
502         m = abs(m);
503 
504         Apfloat inverseModulus = ApfloatMath.inverseRoot(m, 1, m.scale() + Apfloat.EXTRA_PRECISION);
505         a = a.mod(m);
506 
507         if (b.signum() < 0)
508         {
509             // Calculate modular inverse first
510 
            a = modInverse(am);
511             b = b.negate();
512         }
513 
514         Apint two = new Apint(2, b.radix());    // Sub-optimal; the divisor could be some power of two
515 
        Apint[] qr;
516 
517         while ((qr = div(btwo))[1].signum() == 0)
518         {
519             a = modMultiply(aaminverseModulus);
520             b = qr[0];
521         }
522 
523         Apint r = a;
524         qr = div(btwo);
525 
526         while ((b = qr[0]).signum() > 0)
527         {
528             a = modMultiply(aaminverseModulus);
529             qr = div(btwo);
530             if (qr[1].signum() != 0)
531             {
532                 r = modMultiply(raminverseModulus);
533             }
534         }
535 
536         return r;
537     }
538 
539     private static Apint modInverse(Apint aApint m)
540         throws ArithmeticExceptionApfloatRuntimeException
541     {
542         // Extended Euclidean algorithm
543 
        Apint one = new Apint(1, m.radix()),
544               x = Apint.ZERO,
545               y = one,
546               oldX = one,
547               oldY = Apint.ZERO,
548               oldA = a,
549               b = m;
550 
551         while (b.signum() != 0)
552         {
553             Apint q = a.divide(b);
554 
555             Apint tmp = b;
556             b = a.mod(b);
557             a = tmp;
558 
559             tmp = x;
560             x = oldX.subtract(q.multiply(x));
561             oldX = tmp;
562 
563             tmp = y;
564             y = oldY.subtract(q.multiply(y));
565             oldY = tmp;
566         }
567 
568         if (!abs(a).equals(one))
569         {
570             // GCD is not 1
571 
            throw new ArithmeticException("Modular inverse does not exist");
572         }
573 
574         if (oldX.signum() != oldA.signum())
575         {
576             // Adjust by one modulus if sign is wrong
577 
            oldX = oldX.add(copySign(moldA));
578         }
579 
580         return oldX;
581     }
582 
583     /**
584      * 
Factorial function. Uses the default radix.
585      *
586      * 
@param n The number whose factorial is to be calculated. Should be non-negative.
587      *
588      * 
@return <code>n!</code>
589      *
590      * 
@exception java.lang.ArithmeticException If <code>n</code> is negative.
591      * 
@exception java.lang.NumberFormatException If the default radix is not valid.
592      *
593      * 
@since 1.1
594      */

595 
596     public static Apint factorial(long n)
597         throws ArithmeticExceptionNumberFormatExceptionApfloatRuntimeException
598     {
599         return new Apint(ApfloatMath.factorial(nApfloat.INFINITE));
600     }
601 
602     /**
603      * 
Factorial function. Returns a number in the specified radix.
604      *
605      * 
@param n The number whose factorial is to be calculated. Should be non-negative.
606      * 
@param radix The radix to use.
607      *
608      * 
@return <code>n!</code>
609      *
610      * 
@exception java.lang.ArithmeticException If <code>n</code> is negative.
611      * 
@exception java.lang.NumberFormatException If the radix is not valid.
612      *
613      * 
@since 1.1
614      */

615 
616     public static Apint factorial(long nint radix)
617         throws ArithmeticExceptionNumberFormatExceptionApfloatRuntimeException
618     {
619         return new Apint(ApfloatMath.factorial(nApfloat.INFINITEradix));
620     }
621 
622     /**
623      * 
Product of numbers.
624      * 
This method may perform significantly better
625      * 
than simply multiplying the numbers sequentially.<p>
626      *
627      * 
If there are no argumentsthe return value is <code>1</code>.
628      *
629      * 
@param x The argument(s).
630      *
631      * 
@return The product of the given numbers.
632      *
633      * 
@since 1.3
634      */

635 
636     public static Apint product(Apint... x)
637         throws ApfloatRuntimeException
638     {
639         return new Apint(ApfloatMath.product(x));
640     }
641 
642     /**
643      * 
Sum of numbers.<p>
644      *
645      * 
If there are no argumentsthe return value is <code>0</code>.
646      *
647      * 
@param x The argument(s).
648      *
649      * 
@return The sum of the given numbers.
650      *
651      * 
@since 1.3
652      */

653 
654     public static Apint sum(Apint... x)
655         throws ApfloatRuntimeException
656     {
657         return new Apint(ApfloatMath.sum(x));
658     }
659 }
660 
661