GitVita
New repositoryImport repositoryNew issue
Your repositories Your dashboard Sign in

microsoft/qdk

Public

mirrored fromhttps://github.com/microsoft/qdkAvailable

Watch0Fork0Star0
CodeCommitsIssuesPull requestsActionsInsightsSecurity
cd67b36992d2ea20ba330acb56e2e9ac04c11938

Branches

Tags

  • No tags available.
0Branches0Tags
Go to file
Add file
Code

Clone

HTTPS

Download ZIP

qdk/library/std

library/std/math.qs

504lines · modecode

RawDownload
Latest commit unavailable.
1// Copyright (c) Microsoft Corporation.
2// Licensed under the MIT License.
3
4namespace Microsoft.Quantum.Math {
5 open Microsoft.Quantum.Diagnostics;
6 open Microsoft.Quantum.Convert;
7
8 //
9 // Sign, Abs, Min, Max, etc.
10 //
11
12 /// # Summary
13 /// Returns -1, 0 or +1 that indicates the sign of a number.
14 function SignI (a : Int) : Int {
15 if (a < 0) { -1 }
16 elif (a > 0) { +1 }
17 else { 0 }
18 }
19
20 /// # Summary
21 /// Returns -1, 0 or +1 that indicates the sign of a number.
22 function SignD (a : Double) : Int {
23 if (a < 0.0) { -1 }
24 elif (a > 0.0) { +1 }
25 else { 0 }
26 }
27
28 /// # Summary
29 /// Returns -1, 0 or +1 that indicates the sign of a number.
30 function SignL (a : BigInt) : Int {
31 if (a < 0L) { -1 }
32 elif (a > 0L) { +1 }
33 else { 0 }
34 }
35
36 /// # Summary
37 /// Returns the absolute value of an integer.
38 function AbsI (a : Int) : Int {
39 a < 0 ? -a | a
40 }
41
42 /// # Summary
43 /// Returns the absolute value of a double-precision floating-point number.
44 function AbsD (a : Double) : Double {
45 a < 0.0 ? -a | a
46 }
47
48 /// # Summary
49 function AbsL (a : BigInt) : BigInt {
50 a < 0L ? -a | a
51 }
52
53 /// # Summary
54 /// Returns the larger of two specified numbers.
55 function MaxI(a : Int, b : Int) : Int {
56 a > b ? a | b
57 }
58
59 /// # Summary
60 /// Returns the larger of two specified numbers.
61 function MaxD(a : Double, b : Double) : Double {
62 a > b ? a | b
63 }
64
65 /// # Summary
66 /// Returns the larger of two specified numbers.
67 function MaxL (a : BigInt, b : BigInt) : BigInt {
68 a > b ? a | b
69 }
70
71 /// # Summary
72 /// Returns the smaller of two specified numbers.
73 function MinI (a : Int, b : Int) : Int {
74 a < b ? a | b
75 }
76
77 /// # Summary
78 /// Returns the smaller of two specified numbers.
79 function MinD (a : Double, b : Double) : Double {
80 a < b ? a | b
81 }
82
83 /// # Summary
84 /// Returns the smaller of two specified numbers.
85 function MinL(a : BigInt, b : BigInt) : BigInt {
86 a < b ? a | b
87 }
88
89 //
90 // Trigonometric functions
91 //
92
93 /// # Summary
94 /// Represents the ratio of the circumference of a circle to its diameter.
95 ///
96 /// # Ouptut
97 /// A double-precision approximation of the the circumference of a circle
98 /// to its diameter, $\pi \approx 3.14159265358979323846$.
99 ///
100 /// # See Also
101 /// - Microsoft.Quantum.Math.E
102 function PI() : Double {
103 3.14159265358979323846
104 }
105
106 /// # Summary
107 /// Returns the natural logarithmic base to double-precision.
108 ///
109 /// # Output
110 /// A double-precision approximation of the natural logarithic base,
111 /// $e \approx 2.7182818284590452354$.
112 ///
113 /// # See Also
114 /// - Microsoft.Quantum.Math.PI
115 function E() : Double {
116 2.7182818284590452354
117 }
118
119 /// # Summary
120 /// Returns the angle whose cosine is the specified number.
121 function ArcCos (x : Double) : Double {
122 body intrinsic;
123 }
124
125 /// # Summary
126 /// Returns the angle whose sine is the specified number.
127 function ArcSin (y : Double) : Double {
128 body intrinsic;
129 }
130
131 /// # Summary
132 /// Returns the angle whose tangent is the specified number.
133 function ArcTan (d : Double) : Double {
134 body intrinsic;
135 }
136
137 /// # Summary
138 /// Returns the angle whose tangent is the quotient of two specified numbers.
139 function ArcTan2 (y : Double, x : Double) : Double {
140 body intrinsic;
141 }
142
143 /// # Summary
144 /// Returns the cosine of the specified angle.
145 function Cos (theta : Double) : Double {
146 body intrinsic;
147 }
148
149 /// # Summary
150 /// Returns the hyperbolic cosine of the specified angle.
151 function Cosh (d : Double) : Double {
152 body intrinsic;
153 }
154
155 /// # Summary
156 /// Returns the sine of the specified angle.
157 function Sin (theta : Double) : Double {
158 body intrinsic;
159 }
160
161 /// # Summary
162 /// Returns the hyperbolic sine of the specified angle.
163 function Sinh (d : Double) : Double {
164 body intrinsic;
165 }
166
167 /// # Summary
168 /// Returns the tangent of the specified angle.
169 function Tan (d : Double) : Double {
170 body intrinsic;
171 }
172
173 /// # Summary
174 /// Returns the hyperbolic tangent of the specified angle.
175 function Tanh (d : Double) : Double {
176 body intrinsic;
177 }
178
179 /// # Summary
180 /// Computes the inverse hyperbolic cosine of a number.
181 function ArcCosh (x : Double) : Double {
182 Log(x + Sqrt(x * x - 1.0))
183 }
184
185 /// # Summary
186 /// Computes the inverse hyperbolic sine of a number.
187 function ArcSinh (x : Double) : Double {
188 Log(x + Sqrt(x * x + 1.0))
189 }
190
191
192 /// # Summary
193 /// Computes the inverse hyperbolic tangent of a number.
194 function ArcTanh (x : Double) : Double {
195 Log((1.0 + x) / (1.0 - x)) * 0.5
196 }
197
198 //
199 // Sqrt, Log, exp, etc.
200 //
201
202 /// # Summary
203 /// Returns the square root of a specified number.
204 function Sqrt(d : Double) : Double {
205 body intrinsic;
206 }
207
208 /// # Summary
209 /// Returns the natural (base _e_) logarithm of a specified number.
210 function Log(input : Double) : Double {
211 body intrinsic;
212 }
213
214 /// # Summary
215 /// Returns the base-10 logarithm of a specified number.
216 function Log10(input : Double) : Double {
217 Log(input) / Log(10.0)
218 }
219
220 /// # Summary
221 /// Computes the base-2 logarithm of a number.
222 function Lg(input : Double) : Double {
223 Log(input) / Log(2.0)
224 }
225
226 //
227 // Truncation and Rounding
228 //
229
230 /// # Summary
231 /// Returns the integral part of a number.
232 /// For example: Truncate(3.7) = 3; Truncate(-3.7) = -3
233 function Truncate(value : Double) : Int {
234 body intrinsic;
235 }
236
237 internal function ExtendedTruncation(value : Double) : (Int, Double, Bool) {
238 let truncated = Truncate(value);
239 return (truncated, Microsoft.Quantum.Convert.IntAsDouble(truncated) - value, value >= 0.0);
240 }
241
242 /// # Summary
243 /// Returns the smallest integer greater than or equal to the specified number.
244 /// For example: Ceiling(3.1) = 4; Ceiling(-3.7) = -3
245 function Ceiling(value : Double) : Int {
246 let (truncated, remainder, isPositive) = ExtendedTruncation(value);
247 if AbsD(remainder) <= 1e-15 {
248 return truncated;
249 } else {
250 return isPositive ? truncated + 1 | truncated;
251 }
252 }
253
254 /// # Summary
255 /// Returns the largest integer less than or equal to the specified number.
256 /// For example: Floor(3.7) = 3; Floor(-3.1) = -4
257 function Floor(value : Double) : Int {
258 let (truncated, remainder, isPositive) = ExtendedTruncation(value);
259 if AbsD(remainder) <= 1e-15 {
260 return truncated;
261 } else {
262 return isPositive ? truncated | truncated - 1;
263 }
264 }
265
266 /// # Summary
267 /// Returns the nearest integer to the specified number.
268 /// For example: Floor(3.7) = 4; Floor(-3.7) = -4
269 function Round(value : Double) : Int {
270 let (truncated, remainder, isPositive) = ExtendedTruncation(value);
271 if AbsD(remainder) <= 1e-15 {
272 return truncated;
273 } else {
274 let abs = AbsD(remainder);
275 return truncated + (abs <= 0.5 ? 0 | (isPositive ? 1 | -1));
276 }
277 }
278
279 //
280 // Modular arithmetic
281 //
282
283 /// # Summary
284 /// Computes the canonical residue of `value` modulo `modulus`.
285 /// The result is always in the range 0..modulus-1 even for negative numbers.
286 function ModulusI(value : Int, modulus : Int) : Int {
287 Fact(modulus > 0, "`modulus` must be positive");
288 let r = value % modulus;
289 return (r < 0) ? (r + modulus) | r;
290 }
291
292 /// # Summary
293 /// Computes the canonical residue of `value` modulo `modulus`.
294 /// The result is always in the range 0..modulus-1 even for negative numbers.
295 function ModulusL(value : BigInt, modulus : BigInt) : BigInt {
296 Fact(modulus > 0L, "`modulus` must be positive");
297 let r = value % modulus;
298 return (r < 0L) ? (r + modulus) | r;
299 }
300
301 /// # Summary
302 /// Returns an integer raised to a given power, with respect to a given
303 /// modulus. I.e. (expBase^power) % modulus.
304 function ExpModI(expBase : Int, power : Int, modulus : Int) : Int {
305 Fact(power >= 0, "`power` must be non-negative");
306 Fact(modulus > 0, "`modulus` must be positive");
307 Fact(expBase > 0, "`expBase` must be positive");
308
309 // shortcut when modulus is 1
310 if modulus == 1 {
311 return 0;
312 }
313
314 mutable res = 1;
315 mutable expPow2mod = expBase % modulus;
316 mutable powerBits = power;
317
318 while powerBits > 0 {
319 if (powerBits &&& 1) != 0 {
320 // if bit pₖ is 1, multiply res by expBase^(2^k) (mod `modulus`)
321 set res = (res * expPow2mod) % modulus;
322 }
323
324 // update value of expBase^(2^k) (mod `modulus`)
325 set expPow2mod = (expPow2mod * expPow2mod) % modulus;
326 set powerBits >>>= 1;
327 }
328
329 return res;
330 }
331
332 /// # Summary
333 /// Returns an integer raised to a given power, with respect to a given
334 /// modulus. I.e. (expBase^power) % modulus.
335 function ExpModL(expBase : BigInt, power : BigInt, modulus : BigInt) : BigInt {
336 Fact(power >= 0L, "`power` must be non-negative");
337 Fact(modulus > 0L, "`modulus` must be positive");
338 Fact(expBase > 0L, "`expBase` must be positive");
339
340 // shortcut when modulus is 1
341 if modulus == 1L {
342 return 0L;
343 }
344
345 mutable res = 1L;
346 mutable expPow2mod = expBase % modulus;
347 mutable powerBits = power;
348
349 while powerBits > 0L {
350 if (powerBits &&& 1L) != 0L {
351 // if bit pₖ is 1, multiply res by expBase^(2ᵏ) (mod `modulus`)
352 set res = (res * expPow2mod) % modulus;
353 }
354
355 // update value of expBase^(2ᵏ) (mod `modulus`)
356 set expPow2mod = (expPow2mod * expPow2mod) % modulus;
357 set powerBits >>>= 1;
358 }
359
360 return res;
361 }
362
363 /// # Summary
364 /// Returns the multiplicative inverse of a modular integer.
365 ///
366 /// # Description
367 /// This will calculate the multiplicative inverse of a
368 /// modular integer $b$ such that $a \cdot b = 1 (\operatorname{mod} \texttt{modulus})$.
369 function InverseModI(a : Int, modulus : Int) : Int {
370 let (u, v) = ExtendedGreatestCommonDivisorI(a, modulus);
371 let gcd = u * a + v * modulus;
372 Fact(gcd == 1, "`a` and `modulus` must be co-prime");
373 return ModulusI(u, modulus);
374 }
375
376 /// # Summary
377 /// Returns the multiplicative inverse of a modular integer.
378 ///
379 /// # Description
380 /// This will calculate the multiplicative inverse of a
381 /// modular integer $b$ such that $a \cdot b = 1 (\operatorname{mod} \texttt{modulus})$.
382 function InverseModL(a : BigInt, modulus : BigInt) : BigInt {
383 let (u, v) = ExtendedGreatestCommonDivisorL(a, modulus);
384 let gcd = u * a + v * modulus;
385 Fact(gcd == 1L, "`a` and `modulus` must be co-prime");
386 return ModulusL(u, modulus);
387 }
388
389 //
390 // GCD, etc.
391 //
392
393 /// # Summary
394 /// Computes the greatest common divisor of two integers.
395 /// Note: GCD is always positive except that GCD(0,0)=0.
396 function GreatestCommonDivisorI(a : Int, b : Int) : Int {
397 mutable aa = AbsI(a);
398 mutable bb = AbsI(b);
399 while bb != 0 {
400 let cc = aa % bb;
401 set aa = bb;
402 set bb = cc;
403 }
404 return aa;
405 }
406
407 /// # Summary
408 /// Computes the greatest common divisor of two integers.
409 /// Note: GCD is always positive except that GCD(0,0)=0.
410 function GreatestCommonDivisorL(a : BigInt, b : BigInt) : BigInt {
411 mutable aa = AbsL(a);
412 mutable bb = AbsL(b);
413 while bb != 0L {
414 let cc = aa % bb;
415 set aa = bb;
416 set bb = cc;
417 }
418 return aa;
419 }
420
421 /// # Summary
422 /// Returns a tuple (u,v) such that u*a+v*b=GCD(a,b)
423 /// Note: GCD is always positive except that GCD(0,0)=0.
424 function ExtendedGreatestCommonDivisorI(a : Int, b : Int) : (Int, Int) {
425 let signA = SignI(a);
426 let signB = SignI(b);
427 mutable (s1, s2) = (1, 0);
428 mutable (t1, t2) = (0, 1);
429 mutable (r1, r2) = (a * signA, b * signB);
430
431 while r2 != 0 {
432 let quotient = r1 / r2;
433 set (r1, r2) = (r2, r1 - quotient * r2);
434 set (s1, s2) = (s2, s1 - quotient * s2);
435 set (t1, t2) = (t2, t1 - quotient * t2);
436 }
437
438 return (s1 * signA, t1 * signB);
439 }
440
441 /// # Summary
442 /// Returns a tuple (u,v) such that u*a+v*b=GCD(a,b)
443 /// Note: GCD is always positive except that GCD(0,0)=0.
444 function ExtendedGreatestCommonDivisorL(a : BigInt, b : BigInt) : (BigInt, BigInt) {
445 let signA = IntAsBigInt(SignL(a));
446 let signB = IntAsBigInt(SignL(b));
447 mutable (s1, s2) = (1L, 0L);
448 mutable (t1, t2) = (0L, 1L);
449 mutable (r1, r2) = (a * signA, b * signB);
450
451 while r2 != 0L {
452 let quotient = r1 / r2;
453 set (r1, r2) = (r2, r1 - quotient * r2);
454 set (s1, s2) = (s2, s1 - quotient * s2);
455 set (t1, t2) = (t2, t1 - quotient * t2);
456 }
457
458 return (s1 * signA, t1 * signB);
459 }
460
461 /// # Summary
462 /// Finds the continued fraction convergent closest to `fraction`
463 /// with the denominator less or equal to `denominatorBound`
464 /// Using process similar to this: https://nrich.maths.org/1397
465 function ContinuedFractionConvergentI(fraction : (Int, Int), denominatorBound : Int): (Int, Int) {
466 Fact(denominatorBound > 0, "Denominator bound must be positive");
467
468 let (a, b) = fraction;
469 let signA = SignI(a);
470 let signB = SignI(b);
471 mutable (s1, s2) = (1, 0);
472 mutable (t1, t2) = (0, 1);
473 mutable (r1, r2) = (a * signA, b * signB);
474
475 while r2 != 0 and AbsI(s2) <= denominatorBound {
476 let quotient = r1 / r2;
477 set (r1, r2) = (r2, r1 - quotient * r2);
478 set (s1, s2) = (s2, s1 - quotient * s2);
479 set (t1, t2) = (t2, t1 - quotient * t2);
480 }
481
482 return (r2 == 0 and AbsI(s2) <= denominatorBound)
483 ? (-t2 * signB, s2 * signA)
484 | (-t1 * signB, s1 * signA);
485 }
486
487 //
488 // Binary, bits, etc.
489 //
490
491 /// # Summary
492 /// For a non-negative integer `a`, returns the number of bits required to represent `a`.
493 /// NOTE: This function returns the smallest n such that a < 2^n.
494 function BitSizeI(a : Int) : Int {
495 Fact(a >= 0, "`a` must be non-negative.");
496 mutable number = a;
497 mutable size = 0;
498 while (number != 0) {
499 set size = size + 1;
500 set number = number >>> 1;
501 }
502 return size;
503 }
504}
505