text:
Article | Open

# Intrusion of the Kuroshio into the South and East China Seas

• Scientific Reports 7, Article number: 7895 (2017)
• doi:10.1038/s41598-017-08206-4
Accepted:
Published online:

## Abstract

The northward-flowing Kuroshio often intrudes westward and modulates the water masses of the South and East China Seas. These intrusions transcend multiple scales in time and space, which we demonstrate here using various independent data sets. There are two hot spots of intrusion, one in the Luzon Strait and the other off northeast Taiwan, which occur synchronously when the upstream Kuroshio weakens during winter. Beyond seasonal time scales, the two intrusions were not synchronous during 1993–2013. While intrusions into the South China Sea echoed the Pacific Decadal Oscillation, the intrusion northeast of Taiwan decreased markedly before 2002 but regularly reached the shelf thereafter. This change was due to the influence of westward impingements of cyclonic eddies from the open ocean on the Kuroshio main stream in place of anticyclonic eddies. During 1993–2001, decreasing cyclonic eddy impingements moved the Kuroshio away from northeast Taiwan, weakening the Kuroshio intrusion onto the East China Sea shelf. Thereafter, enhanced cyclonic eddy impingement during 2002–2013 weakened the Kuroshio transport, moving it closer to the shelf and enhancing its intrusion into the East China Sea.

## Introduction

The Kuroshio transports water and heat poleward from the tropics. Along its course, it occasionally intrudes into the South China Sea (SCS) through the Luzon Strait and into the East China Sea (ECS) shelf off northeast Taiwan (Fig. 1a). As a result, water masses are exchanged, influencing mass, heat, salinity, and nutrient balances between the Pacific Ocean and SCS and ECS.

Both the seasonal and interannual variabilities of the Kuroshio intrusions have been documented from observations1,2,3,4,5,6,7 and models8, 9. For example, drifters were found to enter the SCS through the Luzon Strait between October and December, but not between July and September2. Recent modeling efforts confirmed that intrusive Kuroshio water is capable of reaching farther westward (~114°E) in winter but is confined to regions east of 118°E in summer. The vertical structure of the intrusive water has also been characterized8, and numerical simulations and limited observations have revealed interannual variations. The winter Kuroshio intrusion into the SCS tends to strengthen during El Niño and weaken during La Niña7. Using Absolute Dynamic Topography, Wu5 demonstrated a strong correlation between Kuroshio intrusion into the SCS and the Pacific Decadal Oscillation (PDO) index. In terms of spatial variations, different types and paths of the Kuroshio intrusion into the SCS have also been identified from satellite altimeter data3, 6.

Shipboard Acoustic Doppler Current Profiler (ADCP) observations between 1995 and 1997 indicated that Kuroshio intrusion onto the ECS shelf off northeast Taiwan was dominant in winter, but that the intrusion was relatively inactive in summer as the Kuroshio migrated seaward4. Further studies indicated that a strengthened Kuroshio migrating seaward coincided with reduced westward intrusions, and a weakened Kuroshio moving shoreward occurred with enhanced westward intrusions. This documented behavior of the Kuroshio was also confirmed using a high-resolution ocean model and appears to transcend weather-band, seasonal and interannual time scales9. Studies to date have only investigated the two intrusion sites separately, and the linkages between them have not been investigated previously.

Two areas, delineated as KIS (115–120°E, 15–20°N)5 and KIE (121.9–122.5°E, 25.5–26.4°N)10 in Fig. 1a, are used to calculate the index of variability in the Kuroshio intrusion into the SCS and ECS, respectively. The KIS index is the monthly mean sea level (in cm) across the region. A low sea level in KIS (west of Luzon Island) is closely correlated with the Kuroshio intrusion into the SCS5. The Kuroshio intrusion current in KIE is generally northeastward with eastward (u) and northward (v) components. The monthly averaged speed index ($u 2 + v 2$ (in cm s−1) in KIE (off northeast Taiwan) accurately reflects the Kuroshio intrusion onto the ECS shelf10. Both intrusion indices were used in the investigations described below.

The Kuroshio is driven by Pacific atmospheric circulation. In the past two decades, pronounced changes in atmospheric circulation and surface winds have been observed globally, such as the strengthening of Pacific trade winds11. Wind-driven circulation in the tropical Pacific Ocean responded to the intensified trade winds with intensified currents and cooling in the eastern Pacific. Furthermore, tropical ocean circulation affected the Kuroshio via the North Pacific subtropical gyre, and fluctuations in the Kuroshio affected intrusions into the SCS and ECS. This study focused on decadal variations in the intrusions, which remain poorly understood. The intrusions into the SCS and onto the ECS shelf displayed both a clear trend and no distinct trend at different times during the last two decades. Different forcing mechanisms are responsible for the two intrusions. In this study, we demonstrated both seasonal and interannual variations as well as residual trends in both intrusions. On interannual and decadal time scales, the two intrusion sites were mostly independent of each other during 1993–2013. We propose a differential forcing mechanism and discuss decadal trends in the other atmospheric factors causing the intrusion of the Kuroshio onto the ECS shelf.

## Seasonal variability

Figure 1a and b show the Argos velocity composites during winter (Dec–Feb) and summer (Jun–Aug), respectively. Clear seasonal patterns were seen in the area of the Luzon Strait. Under the winter northeast monsoon, the northward-flowing Kuroshio meanders to a significant extent into the SCS (Fig. 1a). In summer (Fig. 1b), the corresponding meander is markedly reduced and confined to a region east of 120°E. Observational data clearly demonstrate this contrasting winter–summer behavior. Intrusive Pacific water can be found west of the Luzon Strait in winter but is less clearly identifiable in summer2, 7.

North of the Luzon Strait, the Kuroshio weakens along the east coast of Taiwan in winter (Fig. 1a) but strengthens in summer (Fig. 1b). Based on this behavior, a weakened Kuroshio in winter would lead to its on-shelf encroachment and intrusion, and a strengthened Kuroshio in summer would lead to its seaward migration and decreased intrusion. This anticipated seasonal contrast is too fine in scale to be noticeable off northeast Taiwan in Fig. 1a and b. More specifically, the winter intrusion is widespread over the shelf, but the summer Kuroshio main stream stays off the shelf break (200 m isobath). This seasonal pattern also agrees with the velocity composite of 10-year shipboard ADCP observations12.

Figure 1c and d show the altimeter-based geostrophic velocity in winter and summer, respectively. These seasonal flow patterns and intensity were similar to those derived from Argos drifters (Fig. 1a and b), particularly in the area of the Luzon Strait, including the intrusion position and pathway of the Kuroshio. However, the altimeter data indicated a weaker velocity than did the Argos data. The magnitude of the Argos velocity was closer to that measured by both mooring and shipboard ADCPs; the altimeter-based geostrophic velocity was weaker, in agreement with earlier studies in the area10, 13.

Figure 1e shows the time series of the two Kuroshio intrusion indices derived from the proxy areas KIS and KIE during 1993–2013. The two series are correlated on a seasonal time scale. However, the indices in KIS and KIE showed monthly to interannual variations, and the overall correlation coefficient (R) decreased to a modest value of −0.41 with a significance level of >99%. The lower-than-anticipated R arose from variations other than seasonal. Focusing on only the seasonal variations, the index in KIS peaked, whereas that in KIE reached a minimum in summer, with the reverse pattern observed in winter. Thus, the two intrusions occurred more or less synchronously on seasonal time scales; however, the mechanisms seem to differ. An intense winter monsoon enhanced the Kuroshio intrusion into the SCS14. Northward migration of the North Equatorial Current bifurcation during winter resulted in a weakened Kuroshio off Luzon, also enhancing the intrusion into the SCS5, 15, 16. However, the Kuroshio intrusion onto the ECS shelf may be closely correlated with both the zonal surface heat flux gradient and wind stress in the area. The slower Kuroshio off east Taiwan in winter also promoted on-shelf intrusion10.

## Interannual variability and decadal trend

Since the intrusion at the two sites often takes place during winter, we focused on the winter-only indices (derived by averaging over the period from December to February of the following year) for KIS and KIE to investigate the interannual and decadal variabilities. Figure 2a and b show the resultant KIS and KIE indices in winter during 1993–2013. At a glance, the winter intrusion into the SCS was indicated by a low KIS index in Fig. 2a, and it intensified during 1993–1997 and 2001–2006 but weakened during 1998–2000 and 2007–2013. These episodes reflect the phase change of the PDO (red and blue bars in Fig. 2a). The annual mean PDO index (dashed line) was correlated with the winter-only KIS index during 1993–2013, with an R value of −0.75 above the 99% significance level. Positive PDO generally coincided with a low winter KIS index and enhanced winter intrusions into the SCS, whereas negative PDO led to the opposite pattern. However, a weak correlation exists between the PDO index and KIE winter index (R = 0.29, with a statistical significance below the 90% confidence level) in Fig. 2b. Thus, the PDO generally failed to indicate the winter Kuroshio intrusion onto the ECS shelf.

In addition, the last intrinsic mode functions (IMFs) for the KIS and KIE indices were extracted using the ensemble empirical mode decomposition (EEMD)17 to infer the interannual trends from 1993 to 2013. Also known as the Hilbert-Huang transformation, the EEMD decomposes oscillations with time-varying resolution to preserve nonstationarity and nonlinearity more effectively. It is therefore useful for preserving time-varying time scales intrinsically dictated by the prevailing physics at a specific instant. For each mode of IMF, the number of extrema and zero crossings (after removing the mean) must be equal or differ by no more than one. Each IMF contains linear or nonlinear signals from data and is complete, mostly orthogonal, and adaptive17, 18. Starting from either the KIS or KIE index, EEMD successively removed the highest remaining wavenumber IMF from the time series until the last one (the least wavy IMF) emerged. We focused on the tendency of the last IMF shown in Fig. 2c and d. The last IMF of the winter-only KIS index showed a prominent decadal trend with a >11 cm increase from 1993 to 2013 (Fig. 2c). This increase was within the range reported in previous studies19, 20. Low sea level in KIS generally favors winter Kuroshio intrusion into the SCS. As the low sea level increases, the Kuroshio intrusion in winter should weaken accordingly.

In contrast to the trend displayed by the winter KIS, the last IMF of the winter KIE initially decreased during 1993–2001 and then increased during 2002–2013 (Fig. 2d). Thus, the winter Kuroshio flowing into the ECS tended to stay away from the shelf before 2002 but became more intrusive thereafter. Based on previous studies, the Kuroshio intrusion in the ECS should be inversely correlated with the strength of the upstream Kuroshio off east Taiwan10. To verify whether this inverse relationship transcends interannual time scales, we divided 1993–2013 into two periods, 1993–2001 and 2002–2013. Surface geostrophic currents from satellite altimeter data (AVISO) were used to calculate the linear trend of the upstream Kuroshio intensity for the two periods. The decadal trend was positive during 1993–2001 (Fig. 3a) but negative after 2002 (Fig. 3b), similar to the inverse relationship between incoming Kuroshio strength and on-shelf intrusion on interannual time scales.

Additional data sets also revealed similar tendencies. Surface velocities (near a 15 m depth) were taken from independent data sets that included four ocean reanalysis products (HYCOM, JCOPE-2, GODAS, and ORAS4), each derived from a distinct assimilation method or simulation setting. All data sets consistently showed upstream Kuroshio intensification during the first period (1993–2001) but subsequent weakening (during 2002–2013) with varying magnitudes depending on the data set (Fig. 3c and d). Figure 3c and d correspond to the periods of 1993–2001 and 2002–2013, respectively. For each period, each data set produced a linear fit of the upstream Kuroshio over different areas of coverage (122–123°E and 24–25°N for AVISO, JCOPE2, and HYCOM; 123–124°E and 22–24°N for GODAS; 122–123°E and 23–24°N for ORAS4). Because the data sets differed in resolution, each data set yielded different coverages from the different path width of the mean state of the Kuroshio. However, the area-averaged results from each data set led to the same trends.

## Mechanism

For the residual trend, the two intrusion sites were in phase before 2001, but the southern site became progressively more inactive thereafter, whereas the northern site became more active. The distance between the two intrusion sites is less than 6°. Intervening dynamical process are expected to occur in the area between the two sites, the most likely of which is mesoscale eddy activity. In the area of the Subtropical Countercurrent, the westward-propagating eddies originating from the interior Pacific frequently approach the east coast of Taiwan21. Impinging upon the Kuroshio, the mesoscale eddies are capable of modulating the intensity of the Kuroshio. Cyclonic eddies weaken the Kuroshio, whereas anticyclonic eddies strengthen it22. Therefore, intervention by mesoscale eddies may change the intensity and frequency of intrusion northeast of Taiwan.

Figure 4 shows the time series of the annual eddy characteristics of each polarity following Chelton et al.23 in the proxy area of the Subtropical Countercurrent (122–150°E and 18–25°N), indicated by a rectangle in Fig. 4a. The corresponding number of cyclonic (cold) eddies in Fig. 4b decreased during 1993–2001 but increased thereafter, following similar tendencies to the winter KIE index in Fig. 2d. In contrast, the number of anticyclonic (warm) eddies remained generally constant throughout 1993–2011. The increased cyclonic eddy activity from 2002 to 2013 was consistent with the weakening of the Kuroshio and enhancement of its encroachment and intrusion into the ECS.

In addition to the eddy polarity, the residual statistics of eddy strength and the radius of each polarity were also examined. Among these, only the radius of the cyclonic eddies was statistically significant. A monotonic increase in the radius of the cyclonic eddies during 1993–2013 should further weaken the Kuroshio off east Taiwan.

## Atmospheric connection

Several atmospheric factors may also be relevant to the Kuroshio intrusion onto the ECS shelf. Among these, Chao24 suggested that local northerly wind-induced shelf-ward Ekman drift in winter enhanced the Kuroshio intrusion. Oey et al.25 reported that the seasonal change in the Kuroshio appeared to be related to the surface heat flux gradient over the ECS shelf via the “Joint Effect of Baroclinicity and Relief”. A recent study by Wang et al.26 demonstrated that weakening of the Pacific basin wind stress curl (WSC) reduced the Kuroshio transport and promoted intrusion.

Therefore, we examined the residual trends of local northerly winds, the surface heat flux gradient over the ECS shelf, and basin-wide WSC (data not shown). Among these, the changes in the surface heat flux gradient were not statistically significant during 1993–2013. However, the residual trends of both local northerly winds and basin-wide WSC were fairly constant during 1993–2001 but intensified during 2002–2013. This was consistent with the reinforcement of the effect of the post-2001 increase in cyclonic eddy impingement on the Kuroshio and the associated on-shelf intrusion.

## Conclusion

We analyzed the Kuroshio intrusions into the SCS and onto the ECS shelf using various independent data sets. On seasonal time scales, the intrusions at two sites were generally in phase. However, on interannual and decadal time scales, the two intrusion sites displayed no seasonal correlation during 1993–2013. Different forcing mechanisms were responsible for the two intrusions. The intrusion into the SCS was correlated strongly with the PDO. However, the Kuroshio intrusion onto the ECS shelf was associated with mesoscale eddy activity in the nearby region of the Subtropical Countercurrent. During 1993–2001, decreased cyclonic eddy impingement on the Kuroshio moved the Kuroshio off the shelf break and decreased its intrusion onto the ECS. During 2002–2013, re-strengthened cyclonic eddy impingement led to a weakened Kuroshio, decreasing its on-shelf encroachment and intrusion.

## Methods

### Data set

a. Current data. Ocean surface velocities from Argos drifter data provided by Global Drifter Program (GDP, http://www.aoml.noaa.gov) and Archiving, Validation and Interpretation of Satellite Oceanographic data (AVISO version DT-MADT and DT-MSLA, two sat merged of Ssalto/Duacs, http://www.aviso.altimetry.fr) are plotted in Fig. 1. Argos data provide 15 m depth in-situ measurements of ocean current. The data have been converted to monthly mean with 0.5° spatial resolution before analysis. Validated AVISO data provide surface geostrophic velocities, ADT, and sea level anomaly that are sampled daily with a spatial resolution of 0.25°.

b. Ocean mesoscale eddy data. A statistical dataset of ocean mesoscale eddies is adopted after Chelton et al.23, which provides parameters such as trajectories, polarity, amplitude and radius.

c. Wind stress and heat flux data are from NECPr1 (National Centers for Environmental Prediction/National Center for Atmospheric Research Reanalysis 1)27.

d. The PDO index28 is downloaded from JISAO (Joint Institute for the Study of the Atmosphere and Ocean, http://research.jisao.washington.edu/pdo/).

### Statistical analyses

The significance test of correlation and trend was performed based on the t-test.

### Eddy kinetic energy

The EKE is based on geostrophic calculation29,

$EKE= 1 2 ( U ′ g 2 + V ′ g 2 )$
(1)
$U ′ g =− g f Δ η ′ Δ y$
(2)
$V ′ g = g f Δ η ′ Δ x$
(3)

where $U ′ g$ and $V ′ g$ are the geostrophic velocities, f is the Coriolis parameter, and $η′$ is sea surface height anomaly from AVISO.

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## References

1. 1.

Caruso, M. J., Gawarkiewicz, G. G. & Beardsley, R. C. Interannual variability of the Kuroshio intrusion in the South China Sea. J. Oceanogr. 62, 559–575 (2006).

2. 2.

Centurioni, L. R., Niiler, P. P. & Lee, D.-K. Observations of Inflow of Philippine Sea Surface Water into the South China Sea through the Luzon Strait. J. Phys. Oceanogr. 34, 113–121 (2004).

3. 3.

Nan, F. et al. Identification of different types of Kuroshio intrusion into the South China Sea. Ocean Dyn. 61, 1291–1304 (2011).

4. 4.

Tang, T. Y., Tai, J. H. & Yang, Y. J. The flow pattern north of Taiwan and the migration of the Kuroshio. Cont. Shelf Res. 20, 349–371 (2000).

5. 5.

Wu, C. R. Interannual modulation of the Pacific Decadal Oscillation (PDO) on the low-latitude western North Pacific. Prog. Oceanogr. 110, 49–58 (2013).

6. 6.

Yuan, D., Han, W. & Hu, D. Surface Kuroshio path in the Luzon Strait area derived from satellite remote sensing data. J. Geophys. Res. Oceans 111 (2006).

7. 7.

Yuan, Y. et al. Variation in the Kuroshio intrusion: Modeling and interpretation of observations collected around the Luzon Strait from July 2009 to March 2011. J. Geophys. Res. Oceans 119, 3447–3463 (2014).

8. 8.

Liu, T. et al. Numerical simulation of the Kuroshio intrusion into the South China Sea by a passive tracer. Acta Oceanolog. Sin. 35, 1 (2016).

9. 9.

Wu, C. R., Lu, H. F. & Chao, S. Y. A numerical study on the formation of upwelling off northeast Taiwan. J. Geophys. Res. Oceans 113 (2008).

10. 10.

Wu, C. R., Hsin, Y. C., Chiang, T. L., Lin, Y. F. & Tsui, I. F. Seasonal and interannual changes of the Kuroshio intrusion onto the East China Sea Shelf. J. Geophys. Res. Oceans 119, 5039–5051 (2014).

11. 11.

England, M. H. et al. Recent intensification of wind-driven circulation in the Pacific and the ongoing warming hiatus. Nature Clim. Change 4, 222–227 (2014).

12. 12.

Liang, W. D., Tang, T. Y., Yang, Y. J., Ko, M. T. & Chuang, W. S. Upper-ocean currents around Taiwan. Deep-Sea Res. II 50, 1085–1105 (2003).

13. 13.

Hsin, Y. C., Qu, T. & Wu, C. R. Intra-seasonal variation of the Kuroshio southeast of Taiwan and its possible forcing mechanism. Ocean Dyn. 60, 1293–1306 (2010).

14. 14.

Wu, C. R. & Hsin, Y. C. The forcing mechanism leading to the Kuroshio intrusion into the South China Sea. J. Geophys. Res. Oceans 117 (2012).

15. 15.

He, Y., Cai, S. & Wang, S. The correlation of the surface circulation between the Western Pacific and the South China Sea from satellite altimetry data. Int. J. Remote Sens. 31, 4757–4778 (2010).

16. 16.

Wu, C. R., Wang, Y. L., Lin, Y. F., Chiang, T. L. & Wu, C. C. Weakening of the Kuroshio intrusion into the South China Sea under the global warming hiatus. IEEE J. Sel. Topics Appl. Earth. Observ. Remote Sens. 9, 5064–5070 (2016).

17. 17.

Wu, Z. & Huang, N. E. Ensemble Empirical Mode Decomposition: a noise-assisted data analysis method. Adv. Adapt. data Anal. 1, 1–41 (2009).

18. 18.

Huang, N. E. et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc. R. Soc. London, Ser. A 454, 903–995 (1998).

19. 19.

Cazenave, A. & Remy, F. Sea level and climate: measurements and causes of changes. Wiley Interdiscip. Rev. Clim. Change 2, 647–662 (2011).

20. 20.

Qiu, B., Chen, S., Wu, L. & Kida, S. Wind- versus Eddy-Forced Regional Sea Level Trends and Variability in the North Pacific Ocean. J. Clim. 28, 1561–1577 (2015).

21. 21.

Hwang, C., Wu, C. R. & Kao, R. TOPEX/Poseidon observations of mesoscale eddies over the Subtropical Countercurrent: Kinematic characteristics of an anticyclonic eddy and a cyclonic eddy. J. Geophys. Res. Oceans 109 (2004).

22. 22.

Hsin, Y. C., Qiu, B., Chiang, T. L. & Wu, C. R. Seasonal to interannual variations in the intensity and central position of the surface Kuroshio east of Taiwan. J. Geophys. Res. Oceans 118, 4305–4316 (2013).

23. 23.

Chelton, D. B., Schlax, M. G. & Samelson, R. M. Global observations of nonlinear mesoscale eddies. Prog. Oceanogr. 91, 167–216 (2011).

24. 24.

Chao, S. Y. Circulation of the East China Sea, a numerical study. J. Oceanogr. 46, 273–295 (1990).

25. 25.

Oey, L. Y., Hsin, Y. C. & Wu, C. R. Why does the Kuroshio northeast of Taiwan shift shelfward in winter? Ocean Dyn. 60, 413–426 (2010).

26. 26.

Wang, Y. L., Wu, C. R. & Chao, S. Y. Warming and weakening trends of the Kuroshio during 1993–2013. Geophys. Res. Lett. 43, 9200–9207 (2016).

27. 27.

Kalnay, E. et al. The NCEP/NCAR 40-Year Reanalysis Project. Bull. Amer. Meteor. Soc. 77, 437–471 (1996).

28. 28.

Mantua, N. J., Hare, S. R., Zhang, Y., Wallace, J. M. & Francis, R. C. A Pacific interdecadal climate oscillation with impacts on salmon production. Bull. Amer. Meteor. Soc. 78, 1069–1079 (1997).

29. 29.

Jia, F., Wu, L. & Qiu, B. Seasonal Modulation of Eddy Kinetic Energy and Its Formation Mechanism in the Southeast Indian Ocean. J. Phys. Oceanogr. 41, 657–665 (2011).

## Acknowledgements

This research was supported by the Ministry of Science and Technology, ROC, under grants MOST 104-2611-M-003-002-MY3 and U.S. Office of Naval Research under grant N000141612259.

## Author information

### Affiliations

1. #### National Taiwan Normal University, Taipei, Taiwan

• Chau-Ron Wu
• , You-Lin Wang
•  & Yong-Fu Lin
2. #### University of Maryland, Cambridge, USA

• Shenn-Yu Chao

### Contributions

C.R.W. conceived and designed the study. All authors contributed to the writing and data interpretation. Data analysis: Y.F.L. and Y.L.W. Contextualization: S.Y.C.

### Competing Interests

The authors declare that they have no competing interests.

### Corresponding author

Correspondence to Chau-Ron Wu.


text:
# Definition:Imperial/Length/Inch

## Definition

The inch is an imperial unit of length.

 $\ds$  $\ds 1$ inch $\ds$ $=$ $\ds 2 \cdotp 54$ centimetres $\ds$ $=$ $\ds 0 \cdotp 0254$ metres

This definition is exact, and is how the imperial unit of length is defined.

## Linguistic Note

The word inch derives ultimately from the same root as ounce, that is, from the Latin word uncia, meaning $\dfrac 1 {12}$ part.

In this context, the inch is of course the $\dfrac 1 {12}$ part of a foot.

From the same root, the uncial letters of a mediaeval manuscript are letters which are an inch high.

Despite the gradual migration to the metric system, the word inch still lives on as a rhetorical flourish for a small distance, for example:

Slowly and steadily they inched forward ...
Give him an inch and he'll take a mile.

format: plain


text:
### Archive

Posts Tagged ‘mixed’

## How to simulate multilevel/longitudinal data

I was recently talking with my friend Rebecca about simulating multilevel data, and she asked me if I would show her some examples. It occurred to me that many of you might also like to see some examples, so I decided to post them to the Stata Blog.

### Introduction

We simulate data all the time at StataCorp and for a variety of reasons.

One reason is that real datasets that include the features we would like are often difficult to find. We prefer to use real datasets in the manual examples, but sometimes that isn’t feasible and so we create simulated datasets.

We also simulate data to check the coverage probabilities of new estimators in Stata. Sometimes the formulae published in books and papers contain typographical errors. Sometimes the asymptotic properties of estimators don’t hold under certain conditions. And every once in a while, we make coding mistakes. We run simulations during development to verify that a 95% confidence interval really is a 95% confidence interval.

Simulated data can also come in handy for presentations, teaching purposes, and calculating statistical power using simulations for complex study designs.

And, simulating data is just plain fun once you get the hang of it.

Some of you will recall Vince Wiggins’s blog entry from 2011 entitled “Multilevel random effects in xtmixed and sem — the long and wide of it” in which he simulated a three-level dataset. I’m going to elaborate on how Vince simulated multilevel data, and then I’ll show you some useful variations. Specifically, I’m going to talk about:

1. How to simulate single-level data
2. How to simulate two- and three-level data
3. How to simulate three-level data with covariates
4. How to simulate longitudinal data with random slopes
5. How to simulate longitudinal data with structured errors

### How to simulate single-level data

Let’s begin by simulating a trivially simple, single-level dataset that has the form

$y_i = 70 + e_i$

We will assume that e is normally distributed with mean zero and variance $$\sigma^2$$.

We’d want to simulate 500 observations, so let’s begin by clearing Stata’s memory and setting the number of observations to 500.

. clear
. set obs 500

Next, let’s create a variable named e that contains pseudorandom normally distributed data with mean zero and standard deviation 5:

. generate e = rnormal(0,5)

The variable e is our error term, so we can create an outcome variable y by typing

. generate y = 70 + e

. list y e in 1/5

+----------------------+
|        y           e |
|----------------------|
1. | 78.83927     8.83927 |
2. | 69.97774   -.0222647 |
3. | 69.80065   -.1993514 |
4. | 68.11398    -1.88602 |
5. | 63.08952   -6.910483 |
+----------------------+

We can fit a linear regression for the variable y to determine whether our parameter estimates are reasonably close to the parameters we specified when we simulated our dataset:

. regress y

Source |       SS       df       MS              Number of obs =     500
-------------+------------------------------           F(  0,   499) =    0.00
Model |           0     0           .           Prob > F      =       .
Residual |  12188.8118   499  24.4264766           R-squared     =  0.0000
Total |  12188.8118   499  24.4264766           Root MSE      =  4.9423

------------------------------------------------------------------------------
y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
_cons |   69.89768    .221027   316.24   0.000     69.46342    70.33194
------------------------------------------------------------------------------

The estimate of _cons is 69.9, which is very close to 70, and the Root MSE of 4.9 is equally close to the error’s standard deviation of 5. The parameter estimates will not be exactly equal to the underlying parameters we specified when we created the data because we introduced randomness with the rnormal() function.

This simple example is just to get us started before we work with multilevel data. For familiarity, let’s fit the same model with the mixed command that we will be using later:

. mixed y, stddev

Mixed-effects ML regression                     Number of obs      =       500

Wald chi2(0)       =         .
Log likelihood = -1507.8857                     Prob > chi2        =         .

------------------------------------------------------------------------------
y |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
_cons |   69.89768   .2208059   316.56   0.000     69.46491    70.33045
------------------------------------------------------------------------------

------------------------------------------------------------------------------
Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
sd(Residual) |    4.93737   .1561334      4.640645    5.253068
------------------------------------------------------------------------------

The output is organized with the parameter estimates for the fixed part in the top table and the estimated standard deviations for the random effects in the bottom table. Just as previously, the estimate of _cons is 69.9, and the estimate of the standard deviation of the residuals is 4.9.

Okay. That really was trivial, wasn’t it? Simulating two- and three-level data is almost as easy.

### How to simulate two- and three-level data

I posted a blog entry last year titled “Multilevel linear models in Stata, part 1: Components of variance“. In that posting, I showed a diagram for a residual of a three-level model.

The equation for the variance-components model I fit had the form

$y_{ijk} = mu + u_i.. + u_{ij.} + e_{ijk}$

This model had three residuals, whereas the one-level model we just fit above had only one.


text:
# Toda's theorem

Toda's theorem is a result in computational complexity theory that was proven by Seinosuke Toda in his paper "PP is as Hard as the Polynomial-Time Hierarchy"[1] and was given the 1998 Gödel Prize.

## Statement

The theorem states that the entire polynomial hierarchy PH is contained in PPP; this implies a closely related statement, that PH is contained in P#P.

## Definitions

#P is the problem of exactly counting the number of solutions to a polynomially-verifiable question (that is, to a question in NP), while loosely speaking, PP is the problem of giving an answer that is correct more than half the time. The class P#P consists of all the problems that can be solved in polynomial time if you have access to instantaneous answers to any counting problem in #P (polynomial time relative to a #P oracle). Thus Toda's theorem implies that for any problem in the polynomial hierarchy there is a deterministic polynomial-time Turing reduction to a counting problem.[2]

An analogous result in the complexity theory over the reals (in the sense of Blum–Shub–Smale real Turing machines) was proved by Saugata Basu and Thierry Zell in 2009[3] and a complex analogue of Toda's theorem was proved by Saugata Basu in 2011.[4]

## Proof

The proof is broken into two parts.

• First, it is established that
${\displaystyle \Sigma ^{P}\cdot {\mathsf {BP}}\cdot \oplus {\mathsf {P}}\subseteq {\mathsf {BP}}\cdot \oplus {\mathsf {P}}}$
The proof uses a variation of Valiant–Vazirani theorem. Because ${\displaystyle {\mathsf {BP}}\cdot \oplus {\mathsf {P}}}$ contains ${\displaystyle {\mathsf {P}}}$ and is closed under complement, it follows by induction that ${\displaystyle {\mathsf {PH}}\subseteq {\mathsf {BP}}\cdot \oplus {\mathsf {P}}}$.
• Second, it is established that
${\displaystyle {\mathsf {P}}\cdot \oplus {\mathsf {P}}\subseteq {\mathsf {P}}^{\sharp P}}$

Together, the two parts imply

${\displaystyle {\mathsf {PH}}\subseteq {\mathsf {BP}}\cdot \oplus {\mathsf {P}}\subseteq {\mathsf {P}}\cdot \oplus {\mathsf {P}}\subseteq {\mathsf {P}}^{\sharp P}}$

## References

1. ^ Toda, Seinosuke (October 1991). "PP is as Hard as the Polynomial-Time Hierarchy". SIAM Journal on Computing. 20 (5): 865–877. doi:10.1137/0220053. ISSN 0097-5397.
2. ^ 1998 Gödel Prize. Seinosuke Toda
3. ^ Saugata Basu and Thierry Zell (2009); Polynomial Hierarchy, Betti Numbers and a Real Analogue of Toda's Theorem, in Foundations of Computational Mathematics
4. ^ Saugata Basu (2011); A Complex Analogue of Toda's Theorem, in Foundations of Computational Mathematics


text:
Groupe d’études et de recherche en analyse des décisions

# Two Laplacians for the Distance Matrix of a Graph

## Mustapha Aouchiche et Pierre Hansen

We introduce a Laplacian and a signless Laplacian for the distance matrix of a connected graph, called the distance Laplacian and distance signless Laplacian, respectively. We show the equivalence between the distance signless Laplacian, distance Laplacian and the distance spectra for the class of transmission regular graphs. There is also an equivalence between the Laplacian spectrum and the distance Laplacian spectrum of any connected graph of diameter 2. Similarities between n, as a distance Laplacian eigenvalue, and the algebraic connectivity are established.

, 15 pages

Ce cahier a été révisé en février 2013


text:
# Methylation Helps a Bacterium Adapt to Differing Concentrations

## Bacterial tumbling frequencies remain constant for different attractant concentrations

In the previous lesson, we explored the signal transduction pathway by which E. coli can change its tumbling frequency in response to a change in the concentration of an attractant. But the reality of cellular environments is that the concentration of an attractant can vary across several orders of magnitude. The cell therefore needs to detect not absolute concentrations of an attractant but rather relative changes.

E. coli detects relative changes in its concentration via adaptation to these changes. If the concentration of attractant remains constant for a period of time, then regardless of the absolute value of the concentration, the cell returns to the same background tumbling frequency. In other words, E. coli demonstrates robustness to the attractant concentration in maintaining its default tumbling behavior.

However, our current model is not able to address this adaptation. If the ligand concentration increases in the model, then phosphorylated CheY will plummet and remain at a low steady state.

In this lesson, we will investigate the biochemical mechanism that E. coli uses to achieve such a robust response to environments with different background concentrations. We will then expand the model we built in the previous lesson to replicate the bacterium’s adaptive response.

## Bacteria remember past concentrations using methylation

Recall that in the absence of an attractant, CheW and CheA readily bind to an MCP, leading to greater autophosphorylation of CheA, which in turn phosphorylates CheY. The greater the concentration of phosphorylated CheY, the more frequently the bacterium tumbles.

Signal transduction is achieved through phosphorylation, but E. coli maintains a “memory” of past environmental concentrations through a chemical process called methylation. In this reaction, a methyl group (-CH3) is added to an organic molecule; the removal of a methyl group is called demethylation.

Every MCP receptor contains four methylation sites, meaning that between zero and four methyl groups can be added to the receptor. On the plasma membrane, many MCPs, CheW, and CheA molecules form an array structure. Methylation reduces the negative charge on the receptors, stabilizing the array and facilitating CheA autophosphorylation. The more sites that are methylated, the higher the autophosphorylation rate of CheA, which means that CheY has a higher phosphorylation rate, and tumbling frequency increases.

Note that we now have two different ways that tumbling frequency can be elevated. First, if the concentration of an attractant is low, then CheW and CheA freely form a complex with the MCP, and the phosphorylation cascade passes phosphoryl groups to CheY, which interacts with the flagella and keeps tumbling frequency high. Second, an increase in MCP methylation can also boost CheA autophosphorylation and lead to an increased tumbling frequency.

Methylation of MCPs is achieved by an additional protein called CheR. When bound to MCPs, CheR methylates ligand-bound MCPs faster12, and so the rate of MCP methylation by CheR is higher if the MCP is bound to a ligand.3. Let’s consider how this fact affects a bacterium’s behavior.

Say that E. coli encounters an increase in attractant concentration. Then the lack of a phosphorylation cascade will mean less phosphorylated CheY, and so the tumbling frequency will decrease. However, if the attractant concentration levels off, then the tumbling frequency will flatten, while CheR starts methylating the MCP. Over time, the rising methylation will increase CheA autophosphorylation, bringing back the phosphorylation cascade and raising tumbling frequency back to default levels.

Just as the phosphorylation of CheY can be reversed, MCP methylation can be undone to prevent methylation from being permanent. In particular, an enzyme called CheB, which like CheY is phosphorylated by CheA, demethylates MCPs (as well as autodephosphorylates). The rate of an MCP’s demethylation is dependent on the extent to which the MCP is methylated. In other words, the rate of MCP methylation is higher when the MCP is in a low methylation state, and the rate of demethylation is faster when the MCP is in a high methylation state.3

The figure below adds CheR and CheB to provide a complete picture of the core pathways influencing chemotaxis. To model these pathways and see how our simulated bacterial system responds to different relative attractant concentrations, we will need to add quite a few molecules and reactions to our current model.

The chemotaxis signal-transduction pathway with methylation included. CheA phosphorylates CheB, which methylates MCPs, while CheR demethylates MCPs. Blue lines denote phosphorylation, grey lines denote dephosphorylation, green arrows denote methylation, and red arrows denote demethlyation. Image modified from Parkinson Lab’s illustrations.

## Combinatorial explosion and the need for rule-based modeling

We would like to add the methylation reactions discussed above to the model that we built in the previous lesson and see if this model can replicate the adaptation behavior of E. coli in the presence of a changing attractant concentration. Before incorporating the adaptation mechanisms in our BioNetGen model, we will first describe the reactions that BioNetGen will need.

We begin with considering the MCP complexes. In the phosphorylation tutorial, we identified two components relevant for reactions involving MCPs: a ligand-binding component l and a phosphorylation component Phos. The adaptation mechanism introduces two additional reactions: methylation of the MCP by CheR, and demethylation of the MCP by CheB.

We also need to include binding and dissociation reactions between the MCP and CheR because under normal conditions, most CheR are bound to MCP complexes.4 We will therefore introduce two additional components to the MCP molecules in addition to their phosphorylation components: r (denoting CheR-binding) and Meth (denoting methylation states). In our simulation, we will use three methylation levels (low, medium, and high) rather than five because these three states are most involved in the chemotaxis response to attractants.5

Imagine for a moment that we were attempting to specify every reaction that could take place in our model. To specify an MCP, we would need to tell the program whether it is bound to a ligand (two possible states), whether it is bound to CheR (two possible states), whether it is phosphorylated (two possible states), and which methylation state it is in (three possible states). Therefore, a given MCP has 2 · 2 · 2 · 3 = 24 total states.

Say that we are simulating the simple reaction of a ligand binding to an MCP, which we originally wrote as T + LTL. We now need this reaction to include 12 of the 24 states, the ones corresponding to the MCP being unbound to the ligand. Our previously simple reaction would become 12 different reactions, one for each possible unbound state of the complex molecule T. And if the situation were just a little more complex, with the ligand molecule L having n possible states, then we would have 12n reactions. Imagine trying to debug a model in which we had accidentally incorporated a typo when transcribing just one of these reactions!

In other words, as our model grows, with multiple different states for each molecule involved in each reaction, the number of reactions we need to represent the system grows very fast; this phenomenon is called combinatorial explosion and means that building realistic models of biochemical systems at scale can be daunting.

A major benefit of using a rule-based modeling language such as the one developed by BioNetGen is that it circumvents combinatorial explosion by consolidating many reactions into a single rule. For example, when modeling ligand-MCP binding, we can summarize the 12 different reactions with the rule “a free ligand molecule binds to an MCP that is not bound to a ligand molecule.” In the BioNetGen language, this rule is represented by the same one-line expression as it was in the previous lesson:

LigandReceptor: L(t) + T(l) <-> L(t!1).T(l!1) k_lr_bind, k_lr_dis

We will not bog down the text with a full specification of all the rules needed to add methylation to our model while avoiding combinatorial explosion. If you’re interested in the details, please follow our tutorial.

## Bacterial tumbling is robust to large sudden changes in attractant concentration

In the figures that follow, we plot the concentration over time of each molecule for different values of l0, the initial concentration of ligand. From what we have learned about E. coli, we should see the concentration of phosphorylated CheY (and therefore the bacterium’s tumbling frequency) drop before returning to its original equilibrium. But will our simulation capture this behavior?

First, we add a relatively small amount of attractant, setting l0 equal to 10,000. The system returns so quickly to an equilibrium in phosphorylated CheY that it is difficult to imagine that the attractant has had any effect on tumbling frequency.

Molecular concentrations (in number of molecules in the cell) over time (in seconds) in a BioNetGen chemotaxis simulation with 10,000 initial attractant ligand particles.

If instead l0 is equal to 100,000, then we obtain the figure below. After an initial drop in the concentration of phosphorylated CheY, it returns to equilibrium after a few minutes.

Molecular concentrations (in number of molecules in the cell) over time (in seconds) in a BioNetGen chemotaxis simulation with 100,000 initial attractant ligand particles.

When we increase l0 by another factor of ten to 1 million, the initial drop is more pronounced, but the system returns just as quickly to equilibrium. Note how much higher the concentration of methylated receptors are in this figure compared to the previous figure; however, there are still a significant concentration of receptors with low methylation, indicating that the system may be able to handle an even larger jolt of attractant.

Molecular concentrations (in number of molecules in the cell) over time (in seconds) in a BioNetGen chemotaxis simulation with one million initial attractant ligand particles.

When we set l0 equal to 10 million, we give the system this bigger jolt. Once again, the model returns to its previous CheY equilibrium after a few minutes.

Molecular concentrations (in number of molecules in the cell) over time (in seconds) in a BioNetGen chemotaxis simulation with ten million initial attractant ligand particles.

Finally, with l0 equal to 100 million, we see what we might expect: the steepest drop in phosphorylated CheY yet, but a system that is able to return to equilibrium.

Molecular concentrations (in number of molecules in the cell) over time (in seconds) in a BioNetGen chemotaxis simulation with 100 million initial attractant ligand particles.

Our model, which is built on real reaction rate parameters, provides compelling evidence that the E. coli chemotaxis system is robust to changes in its environment across several orders of magnitude of attractant concentration. This robustness has been observed in real bacteria67, as well as replicated by other computational simulations8.

Aren’t bacteria magnificent?

However, our work is not done. We have simulated E. coli adapting to a single sudden change in its environment, but life is about responding to continual change. So in the next lesson, we will further examine how our simulated bacterium responds in an environment in which the ligand concentration is changing constantly.

If you find chemotaxis biology as interesting as we do, then we suggest the following resources.

1. Amin DN, Hazelbauer GL. 2010. Chemoreceptors in signaling complexes: shifted conformation and asymmetric coupling. Available online

2. Terwilliger TC, Wang JY, Koshland DE. 1986. Kinetics of Receptor Modification - the multiply methlated aspartate receptors involved in bacterial chemotaxis. The Journal of Biolobical Chemistry. Available online

3. Spiro PA, Parkinson JS, and Othmer H. 1997. A model of excitation and adaptation in bacterial chemotaxis. Biochemistry 94:7263-7268. Available online 2

4. Lupas A., and Stock J. 1989. Phosphorylation of an N-terminal regulatory domain activates the CheB methylesterase in bacterial chemotaxis. J Bio Chem 264(29):17337-42. Available online

5. Boyd A., and Simon MI. 1980. Multiple electrophoretic forms of methyl-accepting chemotaxis proteins generated by stimulus-elicited methylation in Escherichia coli. Journal of Bacteriology 143(2):809-815. Available online

6. Shimizu TS, Delalez N, Pichler K, and Berg HC. 2005. Monitoring bacterial chemotaxis by using bioluminescence resonance energy transfer: absence of feedback from the flagellar motors. PNAS. Available online

7. Krembel A., Colin R., Sourijik V. 2015. Importance of multiple methylation sites in Escherichia coli chemotaxis. Available online

8. Bray D, Bourret RB, Simon MI. 1993. Computer simulation of phosphorylation cascade controlling bacterial chemotaxis. Molecular Biology of the Cell. Available online


text:
# My Weblog

## SPOJ 6044. Minimum Diameter Circle [ Enclosing Circle Problem ]

SPOJ 6044. Minimum Diameter Circle problem is also know as enclosing circle problem. Accoriding to Wikipedia ” The smallest circle problem or minimum covering circle problem is a mathematical problem of computing the smallest circle that contains all of a given set of points in the Euclidean plane. The corresponding problem in n-dimensional space, the smallest bounding sphere problem, is to compute the smallest n-sphere that contains all of a given set of points. The smallest circle problem was initially proposed by the English mathematician James Joseph Sylvester in 1857″. Couple of algorithm explained here and implementation Pr. Chrystal algorithm. Accepted Haskell code. You can run this code on ideone.

This algorithm is due to Pr. Chrystal (1885).

The Algorithm

Foremost, observe that the minimal enclosing circle, MEC, is entirely determined by the Convex Hull of a given set of point. The reason is the points of the set touched by the minimal enclosing circle are always on the convex hull of the set.

First step of the algorithm is to compute the convex hull, H, of the points in linear time. Clearly, we can do this since points are kept ordered by x-coordinate. Call the convex hull H and the number of convex hull vertices h. Next step is to pick any side, say S, of the convex hull. The main loop of the algorithm is as follows:

Main Loop

For each vertex of H other than those of S, we compute the angle subtended by S. The minimum such angle occurs at vertex v, and the value of the minimum angle is α (alpha).
1 . IF  α  ≥  90 degrees THEN We are done!
/* The minimal enclosing circle is determined  by the  diametric circle of S */
IF α < 90 degrees THEN
Goto step 2 below.
2.  Since α < 90 degrees, check the remaining vertices of the triangle formed by S and v.
IF no vertices are obtuse THEN We are done!
/* The MEC is determined by the vertices of S and the vertex v */
3. If one of the other angles of the triangle formed by S and v is obtuse, then set S to be the side opposite the obtuse angle and repeat the main loop of the algorithm (i.e., Goto step 1 above).

import Data.List
import qualified Data.Sequence as DS
import Text.Printf
import Data.Function ( on )

data Point a = P a a deriving ( Show , Ord , Eq )
data Turn = S | L | R deriving ( Show , Eq , Ord , Enum  ) -- straight left right

--start of convex hull  http://en.wikipedia.org/wiki/Graham_scan
{--
compPoint :: ( Num  a , Ord a ) => Point a -> Point a -> Ordering
compPoint ( P x1 y1 ) ( P x2 y2 )
| compare x1 x2 == EQ = compare y1 y2
| otherwise = compare x1 x2

findMinx :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ]
findMinx xs = sortBy ( \x  y  -> compPoint  x y  ) xs

compAngle ::(Num a , Ord a ) => Point a -> Point a -> Point a -> Ordering
compAngle ( P x1 y1 ) ( P x2 y2 ) ( P x0 y0 ) = compare ( (  y1 - y0 ) * ( x2 - x0 )  ) ( ( y2 - y0) * ( x1 - x0 ) )

sortByangle :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ]
sortByangle (z:xs) = z : sortBy ( \x y -> compAngle x y z ) xs

convexHull ::( Num a , Ord a )	=> [ Point a ] -> [ Point a ]
convexHull xs = reverse . findHull [y,x]  $ys where (x:y:ys) = sortByangle.findMinx$ xs

findTurn :: ( Num a , Ord a , Eq a ) => Point a -> Point a -> Point a -> Turn
findTurn ( P x0 y0 ) ( P x1 y1 ) ( P x2 y2 )
| ( y1 - y0 ) * ( x2- x0 ) < ( y2 - y0 ) * ( x1 - x0 ) = L
| ( y1 - y0 ) * ( x2- x0 ) == ( y2 - y0 ) * ( x1 - x0 ) = S
| otherwise = R

findHull :: ( Num a , Ord a  )  => [ Point a ] ->   [ Point a ] -> [ Point a ]
findHull [x]  ( z : ys )  = findHull [ z , x ]  ys  --incase of second point  on line from x to z
findHull xs  [] = xs
findHull ( y : x : xs )  ( z:ys )
| findTurn x y z == R = findHull (  x : xs )   ( z:ys )
| findTurn x y z == S = findHull (  x : xs )   ( z:ys )
| otherwise = findHull ( z : y : x : xs  )   ys
--}

--start of monotone hull give .04 sec lead over graham scan

compPoint :: ( Num  a , Ord a ) => Point a -> Point a -> Ordering
compPoint ( P x1 y1 ) ( P x2 y2 )
| compare x1 x2 == EQ = compare y1 y2
| otherwise = compare x1 x2

sortPoint :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ]
sortPoint xs = sortBy ( \ x y -> compPoint x y ) xs

findTurn :: ( Num a , Ord a , Eq a ) => Point a -> Point a -> Point a -> Turn
findTurn ( P x0 y0 ) ( P x1 y1 ) ( P x2 y2 )
| ( y1 - y0 ) * ( x2- x0 ) < ( y2 - y0 ) * ( x1 - x0 ) = L
| ( y1 - y0 ) * ( x2- x0 ) == ( y2 - y0 ) * ( x1 - x0 ) = S
| otherwise = R

hullComputation :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ] -> [ Point a ]
hullComputation [x] ( z:ys ) = hullComputation [z,x] ys
hullComputation xs [] = xs
hullComputation  ( y : x : xs ) ( z : ys )
|  findTurn x y z == R = hullComputation ( x:xs ) ( z : ys )
|  findTurn x y z == S = hullComputation ( x:xs ) ( z : ys )
|  otherwise = hullComputation ( z : y : x : xs ) ys

convexHull :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ]
convexHull xs = final  where
txs = sortPoint xs
( x : y : ys  ) = txs
lhull = hullComputation [y,x] ys
( x': y' : xs' ) = reverse txs
uhull = hullComputation [ y' , x' ] xs'
final = nub $( reverse lhull ) ++ uhull --end of convex hull --start of finding point algorithm http://www.personal.kent.edu/~rmuhamma/Compgeometry/MyCG/CG-Applets/Center/centercli.htm Applet’s Algorithm findAngle :: ( Num a , Ord a , Floating a ) => Point a -> Point a -> Point a -> ( Point a , Point a , Point a , a ) findAngle u@(P x0 y0 ) v@(P x1 y1 ) t@(P x2 y2) | u == t || v == t = ( u , v , t , 10 * pi ) -- two points are same so set the angle more than pi | otherwise = ( u , v, t , ang ) where ang = acos$  ( b + c - a ) / ( 2.0 * ( sqrt $b * c ) ) where sqrDist ( P x y ) ( P x' y' ) = ( x - x' ) ^ 2 + ( y - y' ) ^ 2 [ b , c , a ] = map ( uncurry sqrDist ) [ ( u , t ) , ( v , t ) , ( u , v ) ] findPoints :: ( Num a , Ord a , Floating a ) => Point a -> Point a -> [ Point a ] -> ( Point a , a ) findPoints u v xs | 2 * theta >= pi = circle2Points a b --( a , b , t , theta ) | and [ 2 * alpha <= pi , 2 * beta <= pi ] = circle3Points a b t --( a , b , t , theta ) | otherwise = if 2 * alpha > pi then findPoints v t xs else findPoints u t xs where ( a , b , t , theta ) = minimumBy ( on compare fn ) . map ( findAngle u v )$ xs
fn ( _ , _ , _ , t ) = t
( _ , _ , _ , alpha ) = findAngle v t u  --angle between v u t angle subtended at u by v t
( _ , _ , _ , beta ) = findAngle u t v   -- angle between u v t angle subtended at v by  u t

--end of finding three points algorithm
--find the circle through three points http://paulbourke.net/geometry/circlefrom3/

circle2Points :: ( Num a , Ord a , Floating a ) => Point a -> Point a -> ( Point a , a )
circle2Points ( P x0 y0 ) ( P x1 y1 ) = ( P x y , r ) where
x = ( x0 + x1 ) / 2.0
y = ( y0 + y1 ) / 2.0
r = sqrt $( x0 - x1 ) ^ 2 + ( y0 - y1 ) ^ 2 circle3Points :: ( Num a , Ord a , Floating a ) => Point a -> Point a -> Point a -> ( Point a , a ) --( center , radius ) circle3Points u@(P x1 y1 ) v@(P x2 y2 ) t@(P x3 y3 ) | x2 == x1 = circle3Points u t v | x3 == x2 = circle3Points v u t | otherwise = ( P x y , 2 * r ) where m1 = ( y2 - y1 ) / ( x2 - x1 ) m2 = ( y3 - y2 ) / ( x3 - x2 ) x = ( m1 * m2 * ( y1 - y3 ) + m2 * ( x1 + x2 ) - m1 * ( x2 + x3 ) ) / ( 2 * ( m2 - m1 ) ) y = if y2 /= y1 then ( ( x1 + x2 - 2 * x ) / ( 2 * m1 ) ) + ( ( y1 + y2 ) / 2.0 ) else ( ( x2 + x3 - 2 * x ) / ( 2 * m2 ) ) + ( ( y2 + y3 ) / 2.0 ) r = sqrt$ ( x - x1 ) ^2 + ( y - y1 ) ^ 2

--start of SPOJ code

format::(Num a , Ord a ) => [[a]] -> [Point a]
format xs = map (\[x0 , y0] -> P x0 y0 ) xs

readInt  ::( Num a , Read a ) =>   String -> a

solve :: ( Num a , Ord a , Floating a ) => [ Point a ] -> (  Point a , a )
solve [ P x0 y0 ] = ( P x0 y0 , 0 ) --in case of one point
solve [ P x0 y0 , P x1 y1 ] = circle2Points ( P x0 y0 ) ( P x1 y1 )   -- in case of two points
solve  xs = findPoints x y  t where
t@( x : y : ys )  = convexHull xs

final :: ( Num a , Ord a , Floating a ) => ( Point a , a ) -> a
final ( t , w )
| w == 0 = 0
| otherwise =  w

main = interact $( printf "%.2f\n" :: Double -> String ) . final . solve . convexHull . format . map ( map readInt . words ) . tail . lines  July 24, 2011 ## Monotone Chain Convex Hull Monotone Chain Convex Hull algorithm does not need angular sorting. Wikipedia suggest ” The same basic idea works also if the input is sorted on x-coordinate instead of angle, and the hull is computed in two steps producing the upper and the lower parts of the hull respectively. This modification was devised by A. M. Andrew and is known as Andrew’s Monotone Chain Algorithm. It has the same basic properties as Graham’s Scan but eschews costly comparisons between polar angles “. Also tested this code on SPOJ 3421. Garden Hull and accepted . Also have look at Algorithmist. import Data.List data Point a = P a a deriving ( Show , Ord , Eq ) data Turn = S | L | R deriving ( Show , Ord , Eq , Enum ) --start of monotone convex hull compPoint :: ( Num a , Ord a ) => Point a -> Point a -> Ordering compPoint ( P x1 y1 ) ( P x2 y2 ) | compare x1 x2 == EQ = compare y1 y2 | otherwise = compare x1 x2 sortPoint :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ] sortPoint xs = sortBy ( \ x y -> compPoint x y ) xs findTurn :: ( Num a , Ord a , Eq a ) => Point a -> Point a -> Point a -> Turn findTurn ( P x0 y0 ) ( P x1 y1 ) ( P x2 y2 ) | ( y1 - y0 ) * ( x2- x0 ) < ( y2 - y0 ) * ( x1 - x0 ) = L | ( y1 - y0 ) * ( x2- x0 ) == ( y2 - y0 ) * ( x1 - x0 ) = S | otherwise = R hullComputation :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ] -> [ Point a ] hullComputation [x] ( z:ys ) = hullComputation [z,x] ys hullComputation xs [] = xs hullComputation ( y : x : xs ) ( z : ys ) | findTurn x y z == R = hullComputation ( x:xs ) ( z : ys ) | findTurn x y z == S = hullComputation ( x:xs ) ( z : ys ) | otherwise = hullComputation ( z : y : x : xs ) ys convexHull :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ] convexHull xs = final where txs = sortPoint xs ( x : y : ys ) = txs lhull = hullComputation [y,x] ys ( x': y' : xs' ) = reverse txs uhull = hullComputation [ y' , x' ] xs' final = nub$  ( reverse lhull )  ++  uhull
--end of monotone

--from here on testing part for SPOJ

format::(Num a , Ord a ) => [[a]] -> [Point a]
format xs = map (\[x0 , y0] -> P x0 y0 ) xs

helpSqrt :: (  Floating  a ) => Point a -> Point a ->  a
helpSqrt ( P x0 y0 ) ( P x1 y1 ) =  sqrt  $( x0 - x1 ) ^ 2 + ( y0 - y1 ) ^ 2 solve :: ( Num a , RealFrac a , Floating a ) => [ Point a ] -> String solve xs = ( show . fromIntegral . truncate . snd . foldl ( \( P x0 y0 , s ) ( P x1 y1 ) -> ( P x1 y1 , s + helpSqrt ( P x0 y0 ) ( P x1 y1 ) ) ) ( head xs , 0 )$  ( tail xs  ++ [ head xs ] ) ) ++ "\n"

readInt  ::( Num a , Read a ) =>   String -> a

main = interact $solve . convexHull . format . map ( map readInt . words ) . tail . lines  July 20, 2011 ## Convex Hull in Haskell Implementation Convex hull using Graham scan algorithm and tested it on SPOJ 3421. Garden Hull. Accepted Haskell code.  import Data.List import qualified Data.Sequence as DS data Point a = P a a deriving ( Show , Eq , Ord ) data Turn = S | L | R deriving ( Show , Eq , Ord , Enum ) -- straight left right compPoint :: ( Num a , Ord a ) => Point a -> Point a -> Ordering compPoint ( P x1 y1 ) ( P x2 y2 ) | compare x1 x2 == EQ = compare y1 y2 | otherwise = compare x1 x2 findMinx :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ] findMinx xs = sortBy ( \x y -> compPoint x y ) xs compAngle ::(Num a , Ord a ) => Point a -> Point a -> Point a -> Ordering compAngle ( P x1 y1 ) ( P x2 y2 ) ( P x0 y0 ) = compare ( ( y1 - y0 ) * ( x2 - x0 ) ) ( ( y2 - y0) * ( x1 - x0 ) ) sortByangle :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ] sortByangle (z:xs) = z : sortBy ( \x y -> compAngle x y z ) xs convexHull ::( Num a , Ord a ) => [ Point a ] -> [ Point a ] convexHull xs = reverse . findHull [y,x]$ ys where
(x:y:ys) = sortByangle.findMinx $xs findTurn :: ( Num a , Ord a , Eq a ) => Point a -> Point a -> Point a -> Turn findTurn ( P x0 y0 ) ( P x1 y1 ) ( P x2 y2 ) | ( y1 - y0 ) * ( x2- x0 ) < ( y2 - y0 ) * ( x1 - x0 ) = L | ( y1 - y0 ) * ( x2- x0 ) == ( y2 - y0 ) * ( x1 - x0 ) = S | otherwise = R findHull :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ] -> [ Point a ] findHull [x] ( z : ys ) = findHull [ z , x ] ys --incase of second point on line from x to z findHull xs [] = xs findHull ( y : x : xs ) ( z:ys ) | findTurn x y z == R = findHull ( x : xs ) ( z:ys ) | findTurn x y z == S = findHull ( x : xs ) ( z:ys ) | otherwise = findHull ( z : y : x : xs ) ys --from here on testing part for SPOJ format::(Num a , Ord a ) => [[a]] -> [Point a] format xs = map (\[x0 , y0] -> P x0 y0 ) xs helpSqrt :: ( Floating a ) => Point a -> Point a -> a helpSqrt ( P x0 y0 ) ( P x1 y1 ) = sqrt$  ( x0 - x1 ) ^ 2 + ( y0 - y1 ) ^ 2

solve :: ( Num a , RealFrac a , Floating a  ) => [ Point a ] -> String
solve xs = ( show . fromIntegral . truncate  .  snd . foldl ( \(  P x0 y0  , s )  ( P x1 y1 ) -> ( P x1 y1   , s + helpSqrt  ( P  x0 y0  ) ( P x1 y1 ) ) )  ( head xs , 0 ) $( tail xs ++ [ head xs ] ) ) ++ "\n" readInt ::( Num a , Read a ) => String -> a readInt = read main = interact$ solve . convexHull . format . map  ( map readInt . words ) . tail . lines

July 20, 2011

## Improved Elliptic Curve prime factorisation

Currently i am trying to understand monad. I wrote elliptic curve prime factorisation some times ago and now using monads, i wrote a bit improved elliptic curve prime factorisation . This code factorise the product of prime numbers in range [1000..2000] withing less than second on my slow computer.Algorithm explained is previous post. Modified Haskell code.

--y^2= x^3+ax+b mod n
import Random
import Data.List
import Data.Bits

data Elliptic = Conelliptic Integer Integer Integer deriving ( Eq , Show )
data Point = Conpoint Integer Integer | Identity deriving ( Eq , Show )

powM :: Integer -> Integer -> Integer -> IO Integer
powM a d n = return $powM' a d n where powM' a d n | d == 0 = 1 | d == 1 = mod a n | otherwise = mod q n where p = powM' ( mod ( a^2 ) n ) ( shiftR d 1 ) n q = if (.&.) d 1 == 1 then mod ( a * p ) n else p calSd :: Integer -> IO ( Integer , Integer ) calSd n = return ( s , d ) where s = until ( \x -> testBit ( n - 1) ( fromIntegral x ) ) ( +1 ) 0 d = div ( n - 1 ) ( shiftL 1 ( fromIntegral s ) ) isProbable::Integer->IO Bool isProbable n | n <= 1 = return False | n == 2 = return True | even n = return False | otherwise = calSd n >>= (\( s , d ) -> rabinMiller 0 n s d ) rabinMiller::Integer->Integer->Integer->Integer->IO Bool rabinMiller cnt n s d | cnt>=5 = return True | otherwise = randomRIO ( 2 , n - 2 ) >>= (\a -> powM a d n >>= (\x -> case x == 1 of True -> rabinMiller ( cnt + 1 ) n s d _ -> if any ( == pred n ) . take ( fromIntegral s ) . iterate (\e -> mod ( e^2 ) n )$ x
then rabinMiller ( cnt + 1) n s d
else return False  ) )

addPoints _ Identity p_2 = Left p_2
addPoints _ p_1 Identity = Left p_1
addPoints ( Conelliptic a b n ) ( Conpoint x_p y_p ) ( Conpoint x_q y_q )
| x_p /= x_q = let
[ u , v , d ] = extended_gcd (x_p-x_q) n
s = mod  ( ( y_p - y_q ) * u ) n
x_r = mod ( s * s - x_p - x_q ) n
y_r= mod ( -y_p - s * ( x_r - x_p ) ) n
in case ( ( Conpoint x_r y_r ) , d ) of
( _ , 1 ) -> Left ( Conpoint x_r y_r )
( _ , d' ) -> Right d'
| otherwise = if mod ( y_p + y_q ) n == 0 then Left Identity
else  let
[ u , v , d ] = extended_gcd ( 2*y_p ) n
s = mod  ( ( 3 * x_p * x_p + a ) * u ) n
x_r = mod ( s * s - 2 * x_p ) n
y_r = mod ( -y_p - s * ( x_r - x_p ) ) n
in case ( ( Conpoint x_r y_r ) , d ) of
( _ , 1 )-> Left (Conpoint x_r y_r)
( _ , d' ) -> Right d'

extended_gcd::Integer->Integer->[Integer]
extended_gcd u v= helpGcd u v 0 [ [ 1 , 0 ] , [ 0 , 1 ] ] where
helpGcd u v n m @ ( [ [ a , b ] , [ c , d ] ] )
| v == 0 = if u<0 then [ - ( ( -1 ) ^ n ) * ( m !! 1 !! 1 ) , - ( ( -1 ) ^ ( n + 1 ) ) * ( m !! 0 !! 1 ) , -u ]
else [ ( ( -1 ) ^ n ) * ( m !! 1 !! 1 ) , ( ( -1 ) ^ ( n + 1 ) ) * ( m !! 0 !! 1 ) , u ]
| otherwise = helpGcd u' v' ( n + 1 ) m' where
( q , v' ) = quotRem u v
t = [ [q , 1 ] , [ 1 , 0 ] ]
m' = [ [ q * a + b , a ] , [ q * c + d , c ] ] --mult m t
u' = v

multiEC :: Elliptic -> Point -> Integer -> IO ( Either Point Integer )
multiEC _ _ 0 = return $Left Identity multiEC ecurve point k = return$ helpEC Identity point k
where
helpEC p _ 0 = Left p
helpEC p q n =
case (p',q') of
( Left p'' , Left q'' ) -> helpEC p'' q'' (div n 2)
( Right p'' , _ ) -> Right p''
( _ , Right q'' ) -> Right q''
where
p' = if odd n then addPoints ecurve p q else Left p
q' =  addPoints ecurve q q

dscrmntEC a b = 4 * a * a * a + 27 * b * b

randomCurve :: Integer -> IO [Integer]
randomCurve n =
randomRIO ( 1 , n ) >>=
( \a -> randomRIO ( 1 , n ) >>=
( \u -> randomRIO (1 , n) >>=
(\v -> return [a , mod ( v*v - u*u*u - a*u ) n , n , u , v ] ) ) )

factor :: Integer -> Integer -> IO [Integer]
factor 1 _   = return []
factor n m   =
isProbable n >>= (\x -> if x then  return [n]
else
randomCurve n >>= (\[ a , b , n , u , v ] ->
multiEC ( Conelliptic a b n ) ( Conpoint u v ) m >>=
( \p -> case p of
Left p' -> factor n m
Right p'->  factor ( div n p' ) m  >>= ( \x -> factor p' m  >>= (\y  -> return $x ++ y ) ) ) ) ) solve :: Integer -> IO [ Integer ] solve n = factor n ( foldl lcm 1 [1..10000] ) main = liftM read getLine >>=( \n -> solve n ) >>= print  [user@localhost Haskell]$ time ./Elliptic_Len
14771291995722167329218042370113255785101138530361517229308539087571387624569184860194674818433642463913959873342852255543111135151457209398370030102228858128758106619882182645353368949658870563576166323307466834284670054906895828663228915811739807554988946409285249055824035034910383070200873326731464685909980630202533115332294811841112677175465781631555282945350555352712075708686158031835188673039526232167090966030801645637
[1759,1997,1291,1637,1669,1889,1907,1823,1231,1579,1597,1327,1427,1499,1831,1951,1213,1753,1993,1487,1429,1601,1319,1447,1721,1847,1697,1873,1699,1663,1061,1879,1087,1093,1117,1423,1933,1493,1303,1021,1109,1973,1361,1811,1123,1543,1249,1783,1867,1301,1789,1901,1657,1693,1801,1091,1609,1913,1277,1787,1483,1283,1511,1381,1877,1297,1181,1571,1103,1229,1039,1471,1289,1373,1861,1069,1987,1481,1279,1559,1193,1129,1747,1949,1733,1489,1063,1201,1307,1187,1871,1031,1741,1433,1321,1013,1567,1627,1549,1399,1553,1523,1051,1439,1723,1259,1667,1621,1777,1151,1367,1019,1097,1451,1163,1223,1999,1979,1931,1583,1619,1217,1153,1009,1409,1613,1459,1171,1531,1453,1237,1607,1049,1709,1033]

real	0m16.424s
user	0m1.905s
sys	0m0.009s

July 18, 2011

## SPOJ 6560. N DIV PHI_N (Hard)

SPOJ 6560. N DIV PHI_N (Hard) also accepted using the same idea accept using a bit memoization. Precalculation of product of primes in prime table $<=$ 10^25000 and linear search is enough to get accepted.Its first 6000 primes whose products are sufficient . Accepted Haskell code. I also tried a slightly improved Rabin Miller primality testing with idea give topcoder forum for faster modular exponentiation for large numbers and code with isProbable function got accepted with 3.95 and time limit exceed other time.

import Data.List
import Data.Maybe
import Data.Bits
import qualified Data.ByteString.Char8 as BS

powM::Integer->Integer->Integer->Integer
powM a d n
| d == 0 = 1
| d == 1 = mod a n
| otherwise = mod q n  where
p = powM  ( mod ( a^2 ) n ) ( shiftR d 1 ) n
q = if (.&.) d 1 == 1 then mod ( a * p ) n else p

powL :: Integer -> Integer -> Integer -> Integer
powL a d n
| d <= 10^10 = powM a d n
| otherwise =  mod ( powM ( powL a ( div d 1000000 ) n ) 1000000 n   *  powM a ( mod d 1000000 ) n ) n

isProbable::Integer-> Bool
isProbable n |n<=1 = False
|n==2 = True
|even n = False
|otherwise	= rabinMiller [2,3] n s d where
s = until ( \x -> testBit ( n - 1) ( fromIntegral x ) )  ( +1 ) 0
d = div ( n - 1 ) (  shiftL 1 ( fromIntegral s )  )
--d=until (\x->mod x 2/=0) (div 2) (n-1) --(n-1=2^s*d)
--s=until (\x->d*2^x==pred n) (+1) 0  --counts the power

rabinMiller::[Integer]->Integer->Integer->Integer-> Bool
rabinMiller [] _ _ _ = True
rabinMiller (a:xs) n s d =
if a == n then True
else  case x==1 of
True -> rabinMiller xs n s d
_ -> if any ( == pred n) . take ( fromIntegral s ) . iterate ( \e -> mod (e^2) n ) $x then rabinMiller xs n s d else False where x=powL a d n prime :: [Integer] prime = 2:3:filter isPrime [2*k+1 | k<-[2..]] isPrime :: Integer -> Bool isPrime n = all ((/=0). mod n ).takeWhile ( <= ( truncate.sqrt.fromIntegral$ n) ) $prime anstable ::[Integer] anstable = scanl1 (*) prime helpSolve :: [Integer] -> Integer -> Integer -> Integer helpSolve (x:xs) n ret | x > n = ret | otherwise = helpSolve xs n x solve :: Integer -> BS.ByteString solve 1 = BS.pack "1" solve n = BS.pack.show$ a  where
a = helpSolve ( tail anstable ) n ( head anstable )

main = BS.interact $BS.unlines.map ( solve .readInt ) . tail . BS.lines  July 17, 2011 ## SPOJ 8456. PRIMITIVEROOTS SPOJ 8456. PRIMITIVEROOTS is related to number theory. Product of primitive roots of p ( prime number ) mod p is equal to 1.The product of all primitive roots of prime p $\not=$ 3 is $\equiv 1$ ( mod p ). History of the theory of numbers By Leonard Eugene Dickson on page 183.This problem is not allowed in Haskell. Accepted C++ code. #include<cstdio> #include<iostream> #include<cmath> using namespace std; bool p[10010]; void prime() { for(int i=2 ; i*i <=10000 ; i++ ) if(!p[i]) for(int j = i*i ; j <=10000 ; j += i ) p[j]=1; //for(int i=2; i <= 20 ; i++) if(!p[i]) cout<<i<<" ";cout<<endl; } int main() { prime(); int n , t , x=1 ; for( scanf("%d",&t) ; t>0 ; t--) { scanf("%d",&n); if( n == 3 ) printf("%d:2\n", x++); else if ( !p[n] ) printf("%d:1\n",x++); else printf("%d:NOTPRIME\n",x++); /*if( t != 1 ) printf("\n");*/ } }  July 17, 2011 ## SPOJ Power Sums SPOJ Power Sums accepted after couple of attempts. This problem is similar to SRM 397 SumOfPowers . This discussion on topcoder forum shows couple of ideas to solve this problem using Gaussian elimination , recursive sum and Bernoulli numbers . First i implemented Gaussian elimination mod prime . The idea is $\sum_{i=1}^x i^k = a_{k+1}x^{k+1}+a_kx^k+....+a_1x+a_0$ . For this equation , we have k+2 unknown coefficients so we need k+2 equations to solve this. Putting x from 0 to ( k + 1 ) and we can have k + 2 equation . Solve this system of equation using Gaussian elimination mod prime . $\begin{pmatrix} x_0^{k+1} & x_0^k & . & . & x_0 & 1 \\ x_1^{k+1} & x_1^k & . & . & x_1 & 1 \\ x_2^{k+1} & x_2^k & . & . & x_2 & 1 \\ . & . & . & . & . & . \\ . & . & . & . & . & . \\ x_{k+1}^{k+1} & x_{k+1}^k & . & . & x_{k+1} & 1\end{pmatrix} * \begin{pmatrix} a_{k+1} \\ a_k \\ . \\ . \\ a_1 \\ a_0 \end{pmatrix} = \begin{pmatrix} \sum_{i=1}^{x_0} i^k \\ \sum_{i=1}^{x_1} i^k \\ \sum_{i=1}^{x_2} i^k \\ . \\ .\\ \sum_{i=1}^{x_{k+1}} i^k \end{pmatrix}{\mod p }$ Haskell code for this implementation but it got time limit exceed . This code can also be useb for computation of determinant mod primes and combine the result using chinese remainder theorem. import Data.List extendedGcd :: Integer -> Integer -> ( Integer , Integer ) extendedGcd a b | b == 0 = ( 1 , 0 ) | otherwise = ( t , s - q * t ) where ( q , r ) = quotRem a b ( s , t ) = extendedGcd b r modInv :: Integer -> Integer -> Integer modInv a b | gcd a b /= 1 = error " gcd is not 1 " | otherwise = d where d = until ( > 0 ) ( + b ) . fst.extendedGcd a$ b

adjustwitH::[ Integer ] -> [ Integer ] -> [ Integer ]
adjustwitH ( m:ms ) ( n:ns ) = map ( flip mod 10007 )  $zipWith ( - ) ms$ map ( ( m * modInv n 10007 )* )  ns --m/n mod p modular inverse

echeloN::[[ Integer ]]->[[ Integer ]]
echeloN rs
| null rs || null ( head rs) = rs
| null rs_2 = map (0:) ( echeloN $map tail rs ) --haskell road to logic | otherwise = piv : map ( 0:) ( echeloN rs' ) where rs'= map ( adjustwitH piv ) ( rs_1 ++ rs_3 ) ( rs_1 , rs_2 ) = span ( \( x:xs )-> x == 0 ) rs ( piv:rs_3 ) = rs_2 rows , cols :: [[Integer ]] -> Int rows m = length m cols m = length.head$ m

backS :: [[Integer]] -> [Integer]
backS rs = backS' rs [] where
backS' [] ps = ps
backS' rs ps = backS' rs' ( p:ps ) where
a = ( last rs ) !! ( ( cols rs ) - 2 )
c = ( last rs ) !! ( ( cols rs ) - 1 )
p = mod ( c * modInv a 10007 ) 10007
rs' = eliminate p ( init rs )

eliminate :: Integer -> [[Integer]] -> [[Integer]]
eliminate p rs = map ( simplify p ) rs where
simplify p row = init ( init row ) ++ [ mod ( d - b * p) 10007 ] where
d = last row
b = last ( init row )

helpMod::Integer->Integer->Integer->Integer
helpMod a d n= fun a d where
fun a 0 = 1
fun a 1 = mod a n
fun a 2 = mod (a^2) n
fun a d = let
p=fun (mod (a^2) n) (div d 2)
q=if odd d then mod (a*p) n else p
in mod q n

modS :: Integer -> Integer -> Integer
modS p k = mod ( sum [ helpMod i k 10007 | i <- [1..p]] )  10007

solve :: Integer -> [Integer]
solve k =  backS . echeloN $[ [ helpMod i n 10007 | n <- [ ( k + 1 ) , k .. 0 ] ] ++ [ modS i k ] | i <- [1 .. (k+2)] ] format ::Integer -> [Integer] -> String format k [x] = show x ++ "x^" ++ show k format k t@(x:xs) = ( show x ++ "x^" ++ show k ++ " + ") ++ format ( k - 1 ) xs final :: Integer -> String final k = format ( k + 1) ( solve k ) main = interact$ unlines. map ( final.read ) . lines

My second attempt was recursive formula $( m + 1 )\sum_{k=0}^m = ( n + 1 ) ^{m+1} - \sum_{r=2}^{m+1} ( {{n}\choose{r}} *\sum_{k=0}^n k^{m+1-r} )$ but unfortunately this one also time limit exceed .

import Data.List
import Data.Maybe
import qualified Data.ByteString.Char8 as BS

polySum :: [Integer] -> [Integer] -> [Integer]
polySum ( x : xs ) ys =map ( flip mod 10007 ) $x : zipWith ( + ) xs ys polySub :: [Integer] -> [Integer] -> [Integer] polySub (x:xs) ys = map ( flip mod 10007 )$ x : zipWith ( - ) xs ys

binomM :: Integer -> Integer -> Integer -> Integer
binomM n k p =  pas !! ( fromIntegral n) !! ( fromIntegral k )

pas :: [[Integer]]
pas = iterate ( \x -> map ( flip mod 10007 ) $1 : zipWith (+) x ( tail x ) ++ [1] ) [1] extendedGcd :: Integer -> Integer -> ( Integer , Integer ) extendedGcd a b | b == 0 = ( 1 , 0 ) | otherwise = ( t , s - q * t ) where ( q , r ) = quotRem a b ( s , t ) = extendedGcd b r modInv :: Integer -> Integer -> Integer modInv a b | gcd a b /= 1 = error " gcd is not 1 " | otherwise = d where d = until ( > 0 ) ( + b ) . fst.extendedGcd a$ b

sum' ::[[Integer]] -> [Integer]
sum' xs = foldr ( \x y -> polySum x y ) [0] xs

solve :: Integer -> [[Integer]]
solve m = foldl ( \xs y -> ( polySub ( [ modInv ( y + 1) 10007 * binomM ( y + 1 ) i 10007 | i <- [ 0 .. ( y + 1 ) ] ] ) ( sum' $zipWith map [ ( ( modInv ( y + 1 ) 10007 * binomM ( y + 1 ) r 10007 )* ) | r <- [ 2 .. ( y + 1 ) ] ] xs ) ) : xs ) [[1,1]] [ 1 .. m ]  Finally attempt was using Bernoulli numbers and bingo. Sum is given by $S_m(n) = {1\over{m+1}}\sum_{k=0}^m {m+1\choose{k}} B_k n^{m+1-k}$ and Bernoulli numbers can be computed using $B_m= 1 - \sum_{k=0}^{m-1}\binom mk\frac{B_k}{m-k+1}$.Accepted Haskell code . Learn lot while solving this problem. import Data.List import Data.Maybe import qualified Data.ByteString.Char8 as BS readInt :: BS.ByteString -> Integer readInt = fst . fromJust . BS.readInteger binomM :: Integer -> Integer -> Integer -> Integer binomM n k p = pas !! ( fromIntegral n ) !! ( fromIntegral k ) pas :: [[Integer]] pas = iterate ( \x -> 1 : zipWith (\ a b -> mod ( a + b ) 10007 ) x ( tail x ) ++ [1] ) [1] extendedGcd :: Integer -> Integer -> ( Integer , Integer ) extendedGcd a b | b == 0 = ( 1 , 0 ) | otherwise = ( t , s - q * t ) where ( q , r ) = quotRem a b ( s , t ) = extendedGcd b r modInv :: Integer -> Integer -> Integer modInv a b | gcd a b /= 1 = error " gcd is not 1 " | otherwise = d where d = until ( > 0 ) ( + b ) . fst.extendedGcd a$ b

ber ::  [Integer]
ber  = foldl ( \ xs y ->  if and [ odd y , y /= 1 ] then xs ++ [0] else  xs ++ [ mod ( 1 -  ( sum $zipWith (*) [ binomM y k 10007 * modInv ( y - k + 1 ) 10007 | k <-[ 0 .. pred y ] ] xs ) ) 10007 ] ) [1] [1..300] {--ber is fast than berIt berIt :: [( [Integer] , Integer )] berIt = iterate ( \( xs , y) -> ( map ( flip mod 10007 )$ if and [ odd y , y /= 1 ] then xs ++ [0] else  xs ++ [ 1 -  ( sum $zipWith (*) [ binomM y k 10007 * modInv ( y - k + 1 ) 10007 | k <-[ 0 .. pred y ] ] xs ) ] , y + 1 ) ) ( [1] , 1 ) --} solve :: Integer -> [Integer] solve m = [ mod ( binomM ( m + 1 ) k 10007 * ber !! ( fromIntegral k ) * modInv ( m + 1 ) 10007 ) 10007 | k <- [0..m] ] final :: Integer -> BS.ByteString final k = format ( k + 1 ) ( solve k ) format::Integer -> [Integer] -> BS.ByteString format _ [] = BS.pack "0x^0" format k ( x : xs ) = BS.append ( BS.append ( BS.append ( BS.append ( BS.pack.show$ x ) ( BS.pack  "x^" ) ) ( BS.pack.show $k ) ) ( BS.pack " + " ) ) ( format ( k - 1 ) xs ) main = BS.interact$ BS.unlines . map ( final . readInt )  . BS.lines

Here is more faster way to generate Bernoulli numbers. Replace ber !! ( fromIntegral k ) in solve m with bernoulli k for improved run time.

bernoulli :: Integer -> Integer
bernoulli = genericIndex bernoulli_table
bernoulli_table = map bernoulli' [0..]

bernoulli' :: Integer ->  Integer
bernoulli' 0 = 1
bernoulli' 1 = 1 * modInv 2 10007
bernoulli' n | odd n          = 0
| n mod 6 == 0 = let val = (n+3) * modInv 3 10007 - bernoulliA n (n div 6) in
mod ( val * modInv ( binomM ( n + 3 ) n 10007 ) 10007 ) 10007 --fromIntegral (choose (n+3) n)
| n mod 6 == 2 = let val = (n+3) * modInv 3 10007 - bernoulliA n ((n-2) div 6) in
mod ( val * modInv  ( binomM ( n + 3 ) n 10007 ) 10007 ) 10007 --fromIntegral (choose (n+3) n)
| n mod 6 == 4 = let val = -(n+3) * modInv 6 10007 - bernoulliA n ((n-4) div 6) in
mod (  val * modInv ( binomM ( n + 3 ) n 10007 ) 10007 ) 10007 -- fromIntegral (choose (n+3) n)

bernoulliA :: Integer -> Integer ->  Integer
bernoulliA n m = mod ( sum (map foo [1..m]) ) 10007
where foo j = binomM  (n+3) (n-6*j) 10007  *  bernoulli (n-6*j)

July 15, 2011

## SPOJ 9161. Legendre symbol

SPOJ 9161. Legendre symbol is related to Jacobian symbol.The Jacobi symbol is a generalization of the Legendre symbol. I just submitted my code from Solovay–Strassen primality test and accepted.Accepted Haskell code

import Data.List
import qualified Data.ByteString.Char8 as BS
import Data.Maybe

jacobi::Integral a => a -> a -> a
jacobi a n
| a == 0 = if n == 1 then 1 else 0
| a == 2 && ( mod n 8 == 1 || mod n 8 == 7 ) = 1
| a == 2 && ( mod n 8 == 3 || mod n 8 == 5 ) = -1
| a >= n = jacobi ( mod a n ) n
| even a = jacobi 2 n * jacobi ( div a 2 ) n
| mod a 4 == 3 && mod n 4 == 3 = -jacobi n a
| otherwise = jacobi n a

legend :: [Integer]  -> BS.ByteString
legend (a:p:_) = BS.pack.show $b where b = jacobi a p main = BS.interact$ BS.unlines . map ( legend . map readInt . BS.words ) . tail . BS.lines

July 12, 2011

## SPOJ 7897. Skyline

SPOJ 7897. Skyline is related to Catalan numbers.From wiki ” $C_n$ is the number of permutations of {1, …, n} that avoid the pattern 123 (or any of the other patterns of length 3); that is, the number of permutations with no three-term increasing subsequence”. Accepted Haskell code .

import Data.List
import Data.Ratio
import Data.Maybe
import qualified   Data.ByteString.Char8 as BS
import Data.Array

catlan :: [Integer]
catlan = 1 : zipWith (\c n ->  div ( c * 2 * ( 2 * n + 1 ) )   ( n + 2 )  ) catlan [1..]

catArray_1 ::  Array Int Integer
catArray_1  = listArray ( 0 , 1000 )  catlan

cat :: Array Integer Integer
cat = array (0,1000) ( [( 0 , 1 ) ] ++ [ ( i , div ( cat ! ( i - 1 ) *  ( 4 * i + 2 )  )  ( i + 2 )  ) | i <- [1..1000] ] )

solve :: Integer -> BS.ByteString
solve n = BS.pack.show $mod ( cat ! ( n - 1 ) ) 1000000 main = BS.interact$ BS.unlines. map (  solve.readInt ) . init . BS.lines

July 11, 2011

## SPOJ 1754. Divisor Summation (Hard)

SPOJ 1754. Divisor Summation (Hard) finally accepted in Haskell. I tried this problem in 2007 and almost my 20 solutions in C/C++ and python got time limit exceed. Nothing special , just implemented the pollard rho prime factorisation in Haskell and accepted. Accepted Haskell code.

import Data.List
import Random

helpMod::Integer->Integer->Integer->Integer
helpMod a d n= fun a d where
fun a 1 = mod a n
fun a 2 = mod (a^2) n
fun a d = let
p=fun (mod (a^2) n) (div d 2)
q=if odd d then mod (a*p) n else p
in mod q n

isProbable::Integer-> Bool
isProbable n |n<=1 = False
|n==2 = True
|even n = False
|otherwise	= rabinMiller [2,3,5,7,11,13,17,19] n s d where
d=until (\x->mod x 2/=0) (div 2) (n-1) --(n-1=2^s*d)
s=until (\x->d*2^x==pred n) (+1) 0  --counts the power

rabinMiller::[Integer]->Integer->Integer->Integer-> Bool
rabinMiller [] _ _ _ = True
rabinMiller (a:xs) n s d =  if a == n then True
else  case x==1 of
True -> rabinMiller xs n s d
_ -> if any ( ==pred n) $take (fromIntegral s)$ iterate (\e->mod (e^2) n)  x then rabinMiller xs n s d
else  False
where x=helpMod a d n

fn :: Integer -> Integer -> Integer -> Integer
fn x c n = mod ( x^2 + c ) n

pollard :: Integer -> IO  Integer
pollard n
| even n = return  2
| otherwise =
randomRIO ( 2 , n - 1 ) >>=
(\x' -> randomRIO ( 2 , n - 1 ) >>=
(\c ->  return ( until ( \( _ , _ , d ) -> d /= 1 || d == n )  ( \( x , y , d ) -> ( fn x c n , fn ( fn y c n ) c n , gcd ( fn x c n - fn ( fn y c n ) c n  ) n ) )  ( x' , x' , 1 ) ) >>=
( \( _ , _ , k ) -> return k  ) ) )

factor :: Integer -> IO [ Integer ]
factor 1 = return []
factor n =
(return . isProbable $n ) >>= (\x -> case x of True -> return [n] _ -> pollard n >>= (\d -> case d of 1 -> factor n _ -> factor d >>= (\l -> ( factor.div n$ d ) >>= ( \m -> return $l ++ m ) ) ) ) divSum :: [[Integer]] -> Integer divSum [] = 1 divSum (x:xs) = div ( ( head x )^( genericLength x + 1 ) -1 ) ( head x - 1 ) * divSum xs main = do n <- liftM read getLine let helpfun n cnt | cnt == n = return [] | otherwise = do m <- liftM read getLine lst <- factor m putStrLn.show$ ( divSum.groupBy (\x y -> x == y ).sort \$ lst ) - m
helpfun n (cnt+1)

helpfun n 0

July 8, 2011


text:
1. ## length of curve

find the length of the curve C from t=0 to t pi/2 where x(t)=2sin^3t and y(t)=2cos^3t. set up the correct integral and simplify the integrand.

i think im supposed to use the formula L= sqrt[(dx/dt)^2 + (dy/dt)^2]

for x' i got 36sin^4tcos^2t
y' = -36sin^2tcos^4t

how am i supposed to add(or subtract) them so that it =6sintcost ? that was the answer given for the part inside the integrand. so before rooting im supposed to get 36sin^2tcos^2t but how?

2. Originally Posted by jeph
find the length of the curve C from t=0 to t pi/2 where x(t)=2sin^3t and y(t)=2cos^3t. set up the correct integral and simplify the integrand.

i think im supposed to use the formula L= sqrt[(dx/dt)^2 + (dy/dt)^2]

for x' i got 36sin^4tcos^2t
y' = -36sin^2tcos^4t

how am i supposed to add(or subtract) them so that it =6sintcost ? that was the answer given for the part inside the integrand. so before rooting im supposed to get 36sin^2tcos^2t but how?
$L = \int_0^{\pi /2} dt \, \sqrt{ \left ( \frac{dx}{dt} \right ) ^2 + \left ( \frac{dy}{dt} \right )^2 }$

So
$\frac{dx}{dt} = 6sin^2(3t)cos(3t)$
$\frac{dy}{dt} = -6cos^2(3t)sin(3t)$

Thus
$L = \int_0^{\pi /2} dt \, \sqrt{ \left ( 6sin^2(3t)cos(3t) \right ) ^2 + \left ( -6cos^2(3t)sin(3t) \right ) ^2 }$

$L = \int_0^{\pi /2} dt \, \sqrt{ 36sin^4(3t)cos^2(3t) + 36cos^4(3t)sin^2(3t) }$

$L = 6 \int_0^{\pi /2} dt \, \sqrt{ sin^2(3t)cos^2(3t) (sin^2(3t) + cos^2(3t) ) }$

$L = 6 \int_0^{\pi /2} dt \, sin(3t)cos(3t)$

Can you take it from here?

-Dan

3. i think i might have screwed up on my math typing again...i wanted to say x=2sin^3(t) and y=2cos^3(t). just t not 3t so "sin and cos cubed t" o.o if not then i dont know how to get 3t...

whats the technique to deferentiate those trigs with cubes. can you go
2cos^3(t)
2(cos(t))^3 u=cos(t)
6(u^2)u'
-6cos^2(t)sin(t)

or does the sign switch both times you differentiate? it becomes negative then postive again? or the sign only switches when you differentiate the last cos(t)?

4. Originally Posted by jeph
find the length of the curve C from t=0 to t pi/2 where x(t)=2sin^3t and y(t)=2cos^3t. set up the correct integral and simplify the integrand.

i think im supposed to use the formula L= sqrt[(dx/dt)^2 + (dy/dt)^2]

...
Hello,

$x(t)=2\sin^3(t)$ and $y(t)=2\cos^3(t)$

Then

$\frac{dx}{dt} = 6sin^2(t)cos(t)$ and
$\frac{dy}{dt} = -6cos^2(t)sin(t)$

If you square both sides of these equations you'll get:

$\left(\frac{dx}{dt} \right)^2= 36sin^4(t)cos^2(t)$
$\left( \frac{dy}{dt} \right)^2 = 36cos^4(t)sin^2(t)$

Now calculate the integral:

$L=\int_{0}^{\frac{\pi}{2}} \sqrt{36sin^4(t)cos^2(t) + 36cos^4(t)sin^2(t)} \cdot dt$

$= \int_{0}^{\frac{\pi}{2}} \sqrt{36sin^2(t)cos^2(t)\underbrace{\left( \sin^2(t)+\cos^2(t) \right)}} _{= 1}\cdot dt$

$= 6 \cdot \int_{0}^{\frac{\pi}{2}} |{sin(t)cos(t)} | \cdot dt$

I'll leave the rest for you .

(For confirmation only: I've got L = 3)

5. Originally Posted by jeph
...or the sign only switches when you differentiate the last cos(t)?
Hello,

have a look here:

$x(t) = 2 \sin^3(t)$

$\frac{dx}{dt}= 2 \cdot 3 \cdot \sin^2(t) \cdot(-\cos(t))$

6. Originally Posted by earboth
Hello,

have a look here:

$x(t) = 2 \sin^3(t)$

$\frac{dx}{dt}= 2 \cdot 3 \cdot \sin^2(t) \cdot(-\cos(t))$
shouldnt it be the dy/dt that becomes negative? you differentiate sin(t) right? that just becomes cos(t) not -cos(t)? o.o

7. Originally Posted by jeph
shouldnt it be the dy/dt that becomes negative? you differentiate sin(t) right? that just becomes cos(t) not -cos(t)? o.o
Hello,

of course you are right - but I hope that you caught the method even though I did an ugly mistake. My apologies.

8. Originally Posted by topsquark
$L = \int_0^{\pi /2} dt \, \sqrt{ \left ( \frac{dx}{dt} \right ) ^2 + \left ( \frac{dy}{dt} \right )^2 }$

So
$\frac{dx}{dt} = 6sin^2(3t)cos(3t)$
$\frac{dy}{dt} = -6cos^2(3t)sin(3t)$

Thus
$L = \int_0^{\pi /2} dt \, \sqrt{ \left ( 6sin^2(3t)cos(3t) \right ) ^2 + \left ( -6cos^2(3t)sin(3t) \right ) ^2 }$

$L = \int_0^{\pi /2} dt \, \sqrt{ 36sin^4(3t)cos^2(3t) + 36cos^4(3t)sin^2(3t) }$

$L = 6 \int_0^{\pi /2} dt \, \sqrt{ sin^2(3t)cos^2(3t) (sin^2(3t) + cos^2(3t) ) }$

$L = 6 \int_0^{\pi /2} dt \, sin(3t)cos(3t)$

Can you take it from here?

-Dan
Dan, is there a real reason why you always put $dt$ at the begining of the integral rather than the end? or is it just another one of your personality quirks?

9. Originally Posted by Jhevon
Dan, is there a real reason why you always put $dt$ at the begining of the integral rather than the end? or is it just another one of your personality quirks?
It's a personality quirk. It's also fairly standard in the Physics books I tend to use. There are advantages and disadvantages to its use, but for this level of work I don't see any particular advantage/disadvantage.

-Dan


text:
# Decay rates for elastic-thermoelastic star-shaped networks

Han Z., Zuazua E. Decay rates for elastic-thermoelastic star-shaped networks
Network and Heterogeneous Media, NETW HETEROG MEDIA, Vol. 12, No. 3 (2017), pp. 461-488. DOI: 10.3934/nhm.2017020

Abstract:
This work discusses the asymptotic behaviour of a transmission problem on star-shaped networks of interconnected elastic and thermoelastic rods. Elastic rods are undamped, of conservative nature, while the thermoelastic ones are damped by thermal effects. We analyse the overall decay rate depending of the number of purely elastic components entering on the system and the irrationality properties of its lengths.

First, a sufficient and necessary condition for the strong stability of the thermoelastic-elastic network is given. Then, the uniform exponential decay rate is proved by frequency domain analysis techniques when only one purely elastic undamped rod is present. When the network involves more than one purely elastic undamped rod the lack of exponential decay is proved and nearly sharp polynomial decay rates are deduced under suitable irrationality conditions on the lengths of the rods, based on Diophantine approximation arguments. More general slow decay rates are also derived. Finally, we present some numerical simulations supporting the analytical results.


text:
💬 👋 We’re always here. Join our Discord to connect with other students 24/7, any time, night or day.Join Here!

# Find the absolute maximum and absolute minimum values of $f$ on the given interval.$f(x) = x - 2 \tan^{-1} x$, $[0, 4]$

## Absolute minimum value $1-\dfrac{\pi}{2} \approx -0.5707963268$ which occurs at $x=1$ ; Absolute maximum value $4-2\tan^{-1} \left( 4 \right)$ which occurs at $x=4$

Derivatives

Differentiation

Volume

### Discussion

You must be signed in to discuss.

Lectures

Join Bootcamp

### Video Transcript

let's find the absolute minimum, an absolute maximum values of the function F of x equals x minus two, inverse of the tangent function on x defined on the close interval 04. So first we note that these functions continues over the interval 04. And that because the function universe of the tangent is defined on the real numbers and continue on the real numbers. So it's particularly continue defining continues on the interval 04. Remember that if we uh constrain the tangent function to the domain negative by half by half open interval both sides, then we have in verse well defined. Which is the function that we call tangent to the negative one of X, which is the inverse of the tangent function, which is then like this. So, as we can see, this function is defined and continues under real numbers. Well, so let's say here this is attendant to a negative one of eggs. So f being continuous and defined on a close interval, it attained its extreme values on that interval. And we know that those extreme values can be attained either on the end points of the interval or at critical numbers of the function F. So we get to start by calculating the first derivative of F, Which is 1 -2 times the derivative of the inverse of the tangent function of X is one over one plus six square. So this is equal to 1 -2 over one plus x square. And that's equal to one plus x square as a common denominator. In the numerator we have one plus x square minus two. So the first derivative of f is equal to x square minus one over x square plus one. And uh the denominator we have in this expression is never cereal. In fact it's always positive X square plus one. And so we have that these derivative exists for every uh X in the close interval 04. For that reason we can say that The only critical numbers of f in the interval 04, our total this effects in that interval for which the first derivative is equal to zero. And so we get to solve this equation here and we start with that equation. And because we know this first derivative is equal to the expression here is equivalent to saying that x squared minus one over x square plus one is equal to zero. And this is the same as s square plus one. Sorry, X squared minus one equals zero. That is a square equal one. And from here we know that it's can be negative one, or x can be one, but in this case The Value Native one is not In the interval 04. Then we can say that the only critical number of the function F In the interval 04 is X Equal one. So we have this Google number and the end points of the intervals as candidates for the function attained its extreme values. So we can say that F at the left zero is 0 -2, Inverse of the tangent function at zero. And because the tangent of 00, we get this is equal to zero. And at the writing .4 the function is able to form -2. Universal detergent at four. And is a calculator. We found that these values are approximately equal to one 34 83, 64, 6, 7, 2, 7. And finally F at the critical number x equal one is one minus two. Universe of the tangent at one. And because we know the changing of 5/4 is one. This give us one minus two times by fourth Which is one by half. And that's approximately equal To negative 0.57 0796. Uh 3- 6, 8. Okay, so these are the candidates for string values and functions of the largest of the three values. Is this one here. So it must be The absolute maximum value of the function over the interval 04. And it is occurring at the writing .4. So I will write this here. Absolutely maximum value of the function attained at the right important football. Now, the smallest of the three values is this one here one minus five half. And that attained at the critical number one. So this must be the absolute a minimum value of the function over the interest rate of four. So we write the answer following this observation here. So the the absolute minimum value of the function f uh in sorry on On dangerous 04 is this one here is one minus by health Which is a approximately equal to negative 0.57 0796 3268. And that absolute minimum of course at one. That is a critical number now. 40 absolute maximum Of the function on the interval. 04. He is this one here uh for minus two. Reverse of the tangent function at four which is approximately equal To 1.34836 for 6727. In that absolute maximum value value I forget here while you which of course at de critical number X equal one And I think I hear a road one which is four. It's clear. Look at this, it's just value here. This all here. Is that? Uh so sorry nice one. It was okay. And stopping a bit here. Uh critical number is one. We found it here above. Let's see it again. Critical number is one. So that that value one we have the absolute minimum. So the absolute minimum it's one minus pi half which I presumably equal to this value. And that of course are critical number one was perfect. Now for the absolute maximum value is four. This one here for my institute universal thing that four which presumably equal to that and that occurs. It's not a critical number repeating what it was love. It happens at four. That is the right and point of the interval. So here court at the right endpoint X equal for sorry for that. So this is a great answer and we recall a little bit what we have done here. So the first thing we got to do is to verify that the function is continuous on the clause interval that tells us that the extreme values of that functioning system, that interval are attained either at the end points of the interval or at critical numbers of the function. We calculated the first derivative of the function. We saw that that exists for every argument X in the interval 04. And for that reason the only critical numbers of the function of the in the intervals here for are those were the vax for which The first of the TV- zero we started that equation and there are two solutions native 11 but one of them namely X equal 91 is not in the internal steerforth. So the only critical number of the function f within the interval zero for is X equal one. So we have the critical number one. And the endpoints here. four. So we evaluated at those three points. And with that we found Like the absolute maximum value of the function over the Angels three or 4 which corresponds to the largest value calculated not just of the three values calculated uh Is 4 -2 tangent inverse at four, which is in fact the value that of course at Mexico four. Mhm. So the absolute maximum value Of the function over in two or 3 or four is 4 -2. Universal attended at four that which is approximately equal 1.34836467-7. And it occurs at the writing point of the interview that is at four, That's the largest value we found. And the smallest value, which is uh one minus by half, Which is approximately equal to negative 0.5707963268 Is then the absolute minimum value of the function over the interval 04 and it is attained or of course at the critical number X equal one. And so this is the final answer to the problem.

Derivatives

Differentiation

Volume

Lectures

Join Bootcamp


text:
Applications of Schweitzer's Inequality

Preliminaries

P. Schweitzer proved the following inequality in a 1914 paper. (For a proof, see a separate page.)

For $0\lt m\lt M,\,$ and $x_k\in [m,M],\,$ for $k\in\overline{1,n},$

$\displaystyle \left(\frac{1}{n}\sum_{k=1}^nx_k\right)\left(\frac{1}{n}\sum_{k=1}^n\frac{1}{x_k}\right)\le\frac{(m+M)^2}{4mM}.$

Over the time, several other inequalities have been shown to follow from the above. Tastes may differ, but I found the derivations in Mitrinovic's book both charming and instructive. Three such derivations are reproduced below. I am indebted to N. N. Taleb for bringing Mitrinovic's book to my attention.

Reference

D. S. Mitrinovic, Analytic Inequalities, Springer 1970

Pólya-Szegö Inequality (1925)

Proof

Set $a=0,\,$ $b=\displaystyle \sum_{k=1}^na_kb_k;\,$ set $c_0=0,\,$ $c_t=\displaystyle \sum_{k=1}^ta_kb_k\,$ $t\in\overline{1,n}\,$ and define $f(x)\,$ piecewise: $\displaystyle f(x)=\frac{a_k}{b_k},\,$ $x\in\ (c_{k-1},c_k),\,$ $k\in\overline{1,n}.$

We have $\displaystyle m=\frac{m_1}{M_2}\lt\frac{a_k}{b_k}\lt\frac{M_1}{m_2}=M,\,$ which reduces the required inequality to that of Schweitzer's.

Kantorovich's Inequality (1948)

For $0\lt m\lt\gamma_k\lt M,\,$ and real $u_k,\,$ $k\in\overline{1,n},$

$\displaystyle \left(\sum_{k=1}^{n}\gamma_ku_k^2\right)\left(\sum_{k=1}^{n}\frac{1}{\gamma_k}u_k^2\right)\le\frac{1}{4}\left(\sqrt{\frac{M}{m}}+\sqrt{\frac{m}{M}}\right)^2\left(\sum_{k=1}^{n}u^2_{k}\right)^2.$

Proof

Set $a=0,\,$ $b=\displaystyle \sum_{k=1}^nu_k^2;\,$ set $c_0=0,\,$ $c_t=\displaystyle \sum_{k=1}^tu_k^2\,$ $t\in\overline{1,n}\,$ and define $f(x)\,$ piecewise: $\displaystyle f(x)=\gamma_k,\,$ $x\in\ (c_{k-1},c_k),\,$ $k\in\overline{1,n},\,$ where $0\lt m\lt\gamma_k\lt M,\,$ $k\in\overline{1,n}.$

Now apply Schweitzer's inequality.

Greub-Rheinboldt Inequality (1959)

For $0\lt m_1\lt a_k\lt M_1,\,$ $0\lt m_2\lt b_k\lt M_2,\,$ $k\in\overline{1,n},$

$\displaystyle \left(\sum_{k=1}^{n}a_k^2u_k^2\right)\left(\sum_{k=1}^{n}b_k^2u_k^2\right)\le\frac{(M_1M_2+m_1m_2)^2}{4m_1m_2M_1M_2}\left(\sum_{k=1}^{n}a_kb_ku^2_{k}\right)^2.$

Proof

The Greub-Rheinboldt Inequality is a consequence of that of Kantorovich. This becomes immediately obvious with the substitutions: $\displaystyle \gamma_k=\frac{a_k}{b_k}\,$ and $u_k=\sqrt{a_kb_k}u_k,\,$ where

$0\lt m_1\lt a_k\lt M_1,\,$ $0\lt m_2\lt b_k\lt M_2,\,$ $k\in\overline{1,k}.\,$

We then have

$\displaystyle m=\frac{m_1}{M_2}\le\gamma_k\le\frac{M_1}{m_2}=M$

and Kantorovich's equality yields the Greub-Rheinboldt Inequality.


text:
# DOUBLE SIGMA

Algebra Level 3

the sum ΣΣ(10Cj) (jCi) (where i is less than equal to j which in turn is less then equal to 10. NOTE: i and j are non-negative integers) is

×